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COMPUTER PROGRAM FOR FINITE-DIFFERENCE SOLUTIONS 
OF SHELLS OF REVOLUTION UNDER ASYMMETRIC LOADS 

By Harry G. Schaeffer 
Langley Research Center 

SUMMARY 

A general computer program written in FORTRAN IV symbolic language is pre- 
sented which determines the linear asymmetric bending behavior of an arbitrarily loaded 
elastic thin shell of revolution. The variables are separated by representing the loads, 
displacements, and stresses by Fourier series expansions in the circumferential coor- 
dinate. The resulting set of equations is solved numerically by using finite-difference 
approximations in the meridional direction. A three-layered cross section which is 
symmetric about the middle surface is allowed. The boundary conditions are taken in 
a general form which allows the program to handle elastic restraints by specifying a 
linear combination of edge forces and displacements. The program is described in 
detail and sample calculations are included. 

INTRODUCTION 

The linear behavior of an arbitrary shell of revolution loaded asymmetrically can- 
not, in general, be solved exactly. Numerical techniques for solving the governing 
equations for this general class of problems on high-speed digital computers can, how- 
ever, be employed and several formulations of the numerical solutions of the applicable 
shell equations exist in the literature (see, for example, refs. 1 and 2). The program 
described herein is a mechanization of a finite-difference method substantially as pre- 
sented in reference 1, the analytical formulation being identical. Several features have 
been incorporated in order to make the program as flexible as possible. 

The present report includes sections dealing with analysis, programing techniques, 
and the computer program. Sample problems include the analysis of a cylinder and a 
sphere. In the analysis, the shell material is assumed to be isotropic, and Poisson's 
ratio and the modulus of elasticity are assumed to be constant with respect to the 
directions tangent to the middle surface. The stiffnesses and the Fourier coefficients 
for the surface loads and temperature distribution may have smooth variations in the 
meridional coordinate. The shell may also be of symmetric sandwich construction 



through the thickness provided transverse shear deformations are neglected. Elastic 
constraints at the boundaries are admissible. 

The section describing the computer program is intended to be a user's document 
and contains all the information necessary to prepare input data and subprograms. The 
program is written in FORTRAN IV language for operation in the IBSYS-IBJOB operating 
system (version 13). The program and internal storage requirements utilize approxi- 
mately 25 000 words of memory. 

The program is organized into 25 subprograms to facilitate modifications for 
specific shell shapes and load distributions. A maximum of 500 equal increments along 
the shell meridian is allowed. The program output consists of a problem description 
together with nondimensional displacements, rotations, and moment and force resultants 
in tabular form. 


SYMBOLS 


a 


reference (or characteristic) length 


nondimensional membrane stiffness 


nondimensional bending stiffness 


E 


Young's modulus of elasticity 


reference modulus of elasticity 


e ( n ) e ( n ),e^ Fourier coefficients for membrane strains 

I’ e |fl 

F membrane force 


*4 

f (n) 

u 


Fourier coefficient for transverse shear 

Fourier coefficient for effective (boundary) transverse shear 

shell thickness 


reference thickness 


k^,k^,k^ Fourier coefficients for bending distortion 


2 



L 


last station 


bending moments 

M modified twisting moment (see eq. (12)) 
m^ thermal moment coefficient 

m| n \m^,m^ Fourier coefficients for bending moments 
N total number of shell stations 

N^,N 0 ,N^,Ng^ membrane forces per unit length 

modified membrane shear (see eq. (11)) 
effective (boundary) membrane shear 
n Fourier index 

O first station 

p( n \pj^,p^ Fourier coefficients for loads 

Q^,Q g transverse shear 

effective (boundary) transverse shear 

q,q^q 0 distributed loads in normal, meridional, and circumferential directions, 
respectively 

Rg,Rg principal radii of curvature 

r radial distance from axis of symmetry to shell middle surface 

S total arc length of shell meridian 

s meridional shell coordinate 

T temperature 
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T ( j n) a) 

AT ( 1 n) (?) 

T (n) a,?) 


h (n) 


t (n) . (n) t (n) 
>l 0 ’ l £0 

f(n) 

Ue 

u l ,u e 

u (n) u (n) 

U £ ’ U 0 


w 


w 


(n) 


a 


midplane temperature variation 

temperature gradient per unit thickness normal to the middle surface 
Fourier coefficient for temperature 
cover-plate thickness 
thermal force coefficient 

Fourier coefficients for membrane forces 
Fourier coefficient for effective (boundary) membrane shear 
meridional and circumferential displacements 

Fourier coefficients for meridional and circumferential displacements 
normal displacement 

Fourier coefficient for normal displacement 
coefficient of thermal expansion 



e £ , e fl ,e £0 



0 


meridional increment, 


S 

a(N - 2) 


membrane strains 

coordinate normal to middle surface of shell, positive outward, whose origin 
is at the middle surface 

nondimensional thickness 

circumferential coordinate 
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bending distortions 


K V K e’ K id 



v Poisson's ratio 


5 a 




a rV a 46 1 


J n ) cpW 

t q 


<p 

<P 0 


Subscripts: 


stress components 
reference stress 
rotations 

Fourier components of rotations 

angle between shell center line and normal to shell middle surface 
maximum value of 4> 
nondimensional curvatures 


H horizontal component 

L last station 

max maximum 

O first station 

V vertical component 


Matrices: 


A, B, C, E, F, G ? H, J j A, P 


4X4 matrices 



e,f,g,£,x,y,z 


1x4 matrices 


A prime indicates a derivative with respect to 
FORTRAN names for symbols are as follows: 


Symbol 

FORTRAN NAME 

a 

CHAR 

b 

B 

d 

D 

N 

NMAX 

n 

N 

v 

GAM 

X 

LAM 

V 

NU 

P 

R 

Cl) ^ 
6 

OMT 


OMXI 

A 

DEL 


ANALYTICAL FORMULATION 

The theory and the method of numerical solution are described in detail in refer- 
ence 1. They are summarized in sufficient detail in the text and appendixes so that the 
program can be used without the use of referenced material. 
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Shell Geometry 


The shell geometry and coordinate system for the middle surface of a general 
shell of revolution are shown in figure 1. Any point in the shell may be located by 
specifying the orthogonal coordinates (£, 0,£) where £ = s/a is a nondimensional merid- 
ional coordinate, s is the meridional shell coordinate, a is a reference dimension of 
the shell, 9 is the circumferential coordinate, and £ is a coordinate normal to and 
originating at the middle surface, positive outward. If the shape of the middle surface is 
given by p = p(£) where p = r/a and r is the distance OP, the nondimensional 
principal curvatures can be written as 


a Tl-to^ 1/2 
0 R e 


( 1 ) 


(x) 


= — 
Z " R s 


(y' + y 2 ) 


to . 


where the following expressions have been used: 



( 2 ) 


( 3 ) 


( 4 ) 


y = 


Pl 

p 


( 5 ) 


and a prime denotes differentiation with respect to Finally, the curvatures are 
related by the Codazzi equation 

to 0 ' = y(to£ - to 0 ) (6) 

and the equation 


Pl 1 
P 


-to ..to, 


( 7 ) 
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Equilibrium Equations 


The membrane forces per unit length, transverse force per unit length, moment 
per unit length, and load per unit area are shown in the positive sense in figure 2. The 
displacements and rotations are shown in the positive sense by figure 3; the effective 
boundary force and moment are shown in the positive sense by figure 4. The analysis 
presented in reference 1 is based on the equilibrium equations presented in reference 3. 
The equilibrium equations from reference 1 are written in the following form: 


w, 


4 ( pM l) + U M ee) - P ' M ^| + a [ls( pN {) + U N te) - p ’ N e_ 






+ a^pq£ = 0 


( 8 ) 


JsK) + ilKe) + p,N ?e] + “«^( M e) + ^K») + 

+ 2^("e‘“?)“ie + a2pqe = 0 


0 ) 



l a 

d / \ a / — \ — 

) - p M 0 J 

+ p ~e 

ae( M e) + af( pM ^) + p ' M %e 


( 10 ) 


_ 8 _ 

-ap^N^ + co 0 N ^ + a 2 pq = 0 

where the transverse shear forces and have been eliminated by using the 

moment equilibrium equations. In reference 3 the inplane shear forces and 

as well as the twisting couples and are combined to provide the following 

modified variables: 


N {<. 4 ( N i e + v) + i(i-i)( M 59 - M 9 5 ) 


and 




W -2H« 

which appear in equations (8), (9), and (10). 


Reduction to a Set of Four Ordinary Differential Equations 

It is shown in appendix A and reference 1 that the equilibrium equations (8), (9), 
and (10) may be reduced to a set of four second-order ordinary differential equations. 


(ID 


(12) 
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This set is obtained by expanding the unknowns in appropriate Fourier series and con- 
sidering the strain displacement and stress -strain relations. The resulting set, written 
in matrix form, is given by equation (A28) and is as follows: 


where 


z (n) 


(n) 

(n) 

e 


L m 


(n) 

U 


(13) 


(14) 


and the coefficient matrices E, F, G, and e are defined in appendix A of reference 1. 
These matrices contain the shell stiffnesses and are programed to handle the symmetric 
layered construction shown in figure 5. 


Boundary Conditions 

The set of differential equations which is given by equation (13) is subject to a set 
of boundary conditions at the edges of the shell, or to a set of finiteness conditions at 
points where the radius p becomes zero. The boundary conditions may be any linear 
combination of z and z' at the edge, as is pointed out in reference 1. A set of linear 
relations which is appropriate to the shell problem can be represented by the following 
matrix equation: 

Oy^ n ) + Az^ = l (15) 

where the components of the y^ vector are the Fourier coefficients of the quantities 
N|, Q^, and respectively, and where and are effective membrane 

and transverse shear force intensities, respectively, on the shell boundaries. The ele- 
ments of the £2 and A matrices are specified for a particular boundary condition 
which is to be represented. Examples are given in a later section to illustrate the 
specification of these elements for edge conditions which are representative of those 
which may be found in actual shells. 

The conditions which are to be specified at a shell edge which corresponds to a 
shell closure or a pole (i.e., at points where p becomes zero) are those which are 
presented in references 4 and 5. These conditions are as follows: 
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( 16 ) 




= Ug = w' = m^’ = 0 

(n = 0, 2)" 

+ u g = u^' = w = m^ = 0 

(n = lA 

= u g = w = m^ = 0 

( n = 3) 
J 


COMPUTATIONAL TECHNIQUES 


The set of equations (13) subject to an appropriate set of boundary or finiteness 
conditions is to be solved numerically. The general numerical procedure which is used 
in the present report is very similar to that presented in reference 1. Only the manner 
of handling the boundary conditions is different. The shell meridian showing the finite- 
difference stations used in the present analysis is shown in figure 6. The shell ends O 
and L are halfway between stations 1 and 2 and stations N - 1 and N, respectively. 

The stations associated with points 1 and N lie off the shell and are used to evaluate the 
conditions which exist at the half- stations associated with each end of the shell. In ref- 
erence 1 the need for off-shell points is eliminated by using backward and forward dif- 
ferences at the shell ends; also, half- stations are not used at the ends. The distance 
between stations is equal to A where 

* _ S 
A " a(N - 2) 

and S is the total arc length of the shell. 

At the ends of the shell, vectors Zq, Zq', z^, and z^' are given by the fol- 
lowing relations: 



where the subscripts O and L refer to the first and last shell edges, respectively. 

In the shell interior 2 g i g N - 1, the following central difference approximations 
of the derivatives are used: 
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( 18 ) 


z/ _ ( z i+l ~ z i-l) 




2A 


, „ _ ( z i+l - 2z i + z i-l) 

Ji " A2 




J 


where A is the arc length increment. With these approximations, equation (13) 
becomes a set of linear algebraic equations which can be represented by the following 
matrix form: 


A i z i+1 + B i z i + C i z i-1 = ^ 


(i = 2, 3, . . N - 2, N - 1) 


(19) 


where N is the total number of stations and 


2E< 


A i = F i 

4Ei 

Bi = - + 2AGi 

2Ei 

Ci = ^i-Fi 




g t = 2Ae i 


( 20 ) 


J 


(n) 


It is shown in reference 1 that the vector y v ' which was defined in connection 

(n) * 

with the boundary conditions (see eq. (15)) can be represented in terms of z v ' and 
z^ by the following relation: 


= Hz^' + Jz^ + f 


( 21 ) 


where the nonzero elements of the H, J, and f matrices are given in appendix A of 
reference 1. 

By using equation (21), the boundary conditions can be written entirely in terms of 
z and z\ At the first shell edge, equation (15) becomes 




o( H o z o' + J o z o + f o) + A o z o - 1 o 


( 22 ) 


where the Fourier index has been dropped for convenience. By substituting the applicable 
relations of equations (17) into equation (22) the following relation between z^ and Zg 
is found: 


A 0 Z 2 + B 0 Z 1 = g O 


(23) 
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where 




i! 


B 0 = 


g o = z 

O ■ n o f o v 


( 24 ) 


Because of the character of the system of equations to be solved a general result 
for z^ in terms of can be written as follows: 


V 


- p i z i+l + x i 


(25) 


It would appear to be expedient to determine Pj and directly by solving the matrix 
equation (23) and equating coefficients of equation (25); however, the matrix Bq is 
singular for the case where displacements and rotation are prescribed. The solution is 
thus obtained by solving first for from the equilibrium equation written at station 2 
as follows: 

Z 1 = C 2 ( g l " A 2 Z 3 " B 2 z 3) (26) 

By substituting equation (26) into equation (23), Z 2 is determined as follows: 

z 2 = ( B 0 C 2 1b 2 " A o) (“ B 0 C 2 1 a 2 z 3 + B 0 C 2 lg 2 " g o) (27) 

The matrices Pg and Xg are seen by inspection to be as follows: 

P 2 ' ( B 0 C 2 1b 2 - A o) 1b O C 2 1a 2 
x 2 ' ( B 0 C 2 1b 2 " A o) ( B 0 C 2 lg 2 ' g O 



The general expressions for the Pj and x A matrices, which are valid for 
2 = i = N- 1, are found by substituting equation (25) into equation (19) and are 


p > = ( B i - c i p i-i)' lA i 

x i - ( B i - C i P i-l)" 1 (s i - C i x i-l) 


(29) 


i 

It is seen that, after calculating P 2 and X 2 from equation (28), Pj and x^ can be 
calculated at each interior shell station by successively applying equation (29). 
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In order to calculate the solution vector at each station, must be determined. 

The vector is found by simultaneously solving the boundary equation (15) and equa- 

tion (25) for i = N - 1. The result is as follows: 


Z N = 


n Ll A T/ IT 


£2-r 


^ L " fi L f L 


„ ) + 0k 

a + 2 r~ 


A 2/ 2 


x 


N- 1/ 


N-ll 


(30) 


The solution can now be determined at each station by successive application of equa- 
tion (25). After z^ (i = 1, . . N) have been calculated all the other shell quantities 
of interest can be determined. The quantities which are calculated and presented as 
program output are the nondimensional Fourier coefficients t^, t^, f^, cp^, m 0 , 
m and t 0 . The quantity f^ is the transverse shear and is given by the relation 


nm 


f | = y(m^-m 0 ) + m^ + ? 


if 


(31) 


The remaining quantities are defined in terms of z * in reference 1. 


COMPUTER PROGRAM 


The set of equations and the associated boundary conditions presented in refer- 
ence 1 and described in other sections of this report have been written in FORTRAN IV 
language for operation in an IBSYS-IBJOB operating system (version 13). In deriving 
the set of equations the unknowns are expressed as suitable Fourier expansions; thus, 
the program calculates a set of Fourier coefficients for each value of the Fourier index. 
These coefficients may then be summed by using the appropriate Fourier series to give 
the total value of the variable as a function of the circumferential coordinate. 

In developing any computer program a basic decision must be made regarding its 
organization. On the one hand, many of the important parameters characterizing the 
problem to be solved may be predefined, limiting to a great degree the burden imposed 
on the user of providing input information but at the same time limiting the flexibility of 
the program. On the other hand, the program may be kept general by requiring that the 
user prepare subprograms which define the important parameters as mathematical 
functions. The latter choice has been made herein. The user must prepare compatible 
subprograms which define the geometrical properties of the middle surface and which 
specify the loads, temperature, and stiffnesses along the shell meridian. In addition, 
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input data must be provided for the edge conditions and for the user-prepared subpro- 
grams. A detailed description of input requirements is presented in a later section. 

Program Organization 

The flow diagram for the main program is presented in figure 7. It can be seen 
that the main program causes input data to be read and then calls for a number of sub- 
routines in a logical order. As an aid in reading the flow chart, a list of subroutines 
and their description is presented in table 1. As diagramed in figure 7, the calculations 
and logical operations are carried out sequentially in the following order: 

I. Read input data and program control words 

(1) Calculate geometric parameters of reference surface if the control word 

called IND3 * 0. 

(2) List problem description and input data. 

II. Calculate the P and x matrix and vector, respectively, at the shell sta- 
tion i = 2. This step utilizes subprograms called HFJ, FORCE, ABCG, 
and INIT. 

in. Calculate and store P and x at each shell station through i = N - 1. This 
step utilizes EFG, FORCE, and PANDX. 

IV. Calculate and store the vector z^. This step utilizes HFJ and FINAL. 

V. Calculate and store the vector zj associated with each shell station starting 
at i = N and proceeding to i = 1. This step utilizes EQ73. 

VI. Calculate output quantities t^, t^, t^, f^, m 0 , m <p ^ and print 
results. This step utilizes STRESS AND OUTPUT. 

A complete list of the FORTRAN IV source deck for the main program, SHELLS, 
and all the subroutines is presented in appendix B. A glossary of the FORTRAN names 
of the variables is given in table 2. 

Input Data 

The format for the required and optional input data cards is now given. These 
cards must be in the sequence given. 

1. Shell properties and program control indicators: (FORMAT 5F6.3, 614) 

NU, TKN, CHAR, N, EALSIG, NO, IND1, IND2, IND3, IND4, IND5 
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where 


NU 

TKN 

CHAR 

N 

EALSIG 

NO 

IND1 

= 0 
* 0 

IND2 

IND3 

= 0 
* 0 


Poisson's ratio (floating point) 

reference thickness to be selected by the user and is used internally to 
nondimensionalize the results 

reference shell dimension to be selected by user and is used internally to 
nondimensionalize the results 

Fourier index (floating point) 

Ecu 

— , thermal coefficient used only if thermal stresses are calculated 

u o 

problem number (fixed point) 

an integer value for option of printing geometrical properties of reference 
surface as a part of the output: 

geometrical data are printed 

geometrical data not printed 

an integer indicating the number of words to be read into an array called 
CONST which is described later. The number of words is equal to the 
numerical value of IND2 (IND2 s 100) 

the call by the main program for the calculation of shell geometric param- 
eters is conditional depending on the numerical value of this indicator as 
follows: 

the shell geometric data which have been previously calculated and stored 
are to be used. This allows several load distributions to be applied to the 
same shell configuration without the necessity of recalculating the geomet- 
ric parameters. 

the shell geometric parameters are to be calculated. The value of IND3 must 
be other than zero for the first problem of a sequence of problems and for 
any succeeding problem in a sequence when a change of geometry is 
desired. When IND3 * 0, additional input cards are required as described 
under item 4. 
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IND4 


boundary condition matrix control: 


= 0 indicates that only diagonal elements of £2 and A matrices are to be read. 
The off-diagonal terms are automatically set to zero. 

* 0 indicates that entire £2 and A matrices are to be read. 

IND5 an integer value of 0, 1, 2, or 3: 

= 0 no poles 

= 1 pole at £ = 0 

= 2 poles at £ = 0, £ max 

= 3 pole at | = 

2. Problem description; 72 characters of alphanumeric information in order to 
describe briefly the problem being run. 

3. Variable parametric data: Data card or cards appear only if the indicator 
IND2 #0. If IND2 + 0, HSTD2 numbers are to be provided by using a maximum of 25 cards 
having FORMAT 4E16.8. These data are stored in the CONST array and are available 

in COMMON location BL14. 

4. Geometric data: The following card or cards appear only if the indicator 

IND3 ^0. If IND3 * 0, the following data are to be presented: 

a. FORMAT (214) 

NMAX, FREQ 

where 

NMAX total number of stations (NMAX S 502) 

FREQ integer which controls the frequency for printing numerical 

results. Results are printed at every FREQth station. 

b. Cards are inserted as required for calculation of shell geometric prop- 

erties. The format of such cards is to be specified by the INPUT 
subroutine. 

5. Boundary Condition Cards (FORMAT 4E16.8): Two indicators, IND4 and IND5, 

are used to specify first the need for and second the format of the data which specify the 
boundary condition at each end of the shell. For the case where a pole exists on the 
shell, as indicated by an appropriate value of IND5, the appropriate finiteness conditions 
are generated within the program and no cards are to be supplied for these conditions. 

If, however, the value of IND5 indicates that no pole is present, boundary condition cards 
must be supplied in one of two formats depending on the value of IND4. 
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a. IND4 = 0 


signifies that only the diagonal elements of and A are 
required. The diagonal elements of the boundary condition 
matrices and the vector l are to appear in the following 


order: 




(i) 

^11’ ^22’ 

^33 > 

^44 

(2) 

A ll> a 22> 

a 33> 

a 44 

(3) 

^2’ ^3’ 

1 4 



b. IND4 ± 0 signifies that the entire arrays O and A are required. 

This allows the specification of elastic boundary conditions 
or the specification of a force or displacement in a direc- 
tion other than an intrinsic coordinate direction. The ele- 
ments of these arrays and the l -vector are to appear in the 
following order: 


(1) 

%1> 

n 21> 

n 31> 

n 41 

(2) 

^12> 

^22’ 

fi 32> 

n 42 

(3) 

n 13> 

n 23> 

n 33> 

CO 

cf 

(4) 

^14> 

^24’ 

fi 34> 

n 44 

(5) 

A ll> 

a 21 > 

A 31’ 

A 41 

(6) 

A 12’ 

A 22> 

A 32> 

A 42 

(7) 

A 13’ 

A 23> 

A 33> 

> 

4^ 

CO 

(8) 

a 14> 

a 24> 

> 

CO 

4^ 

A 44 

(9) 

l V l 

2> l 3> 

h 



It should be noted that the cards described must be in the indicated order with those 
cards associated with the edge £ = 0 preceding those associated with the edge 
£ = In addition, if IND4 * 0, the full arrays O and A are required at each 

edge unless, of course, the need of either or both sets of cards is negated by an appro- 
priate value of IND5. If IND5 = 2, indicating a pole at each end of the shell, no boundary 
condition cards are to be supplied; all the necessary conditions are generated by the 
program. 

6. Punch control card (FORMAT 1114): The option of punching any or all of the 
11 output variables is controlled by this card. If it is desired to punch data, a one is 
punched in the appropriate column of the punch control card as follows: 
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Variable 


Column 


H 

4 

t e 

8 

t ^e 

12 

h 

16 


20 

m # 

24 


28 


32 

u e 

36 

w 

40 


44 


If no punching is desired, a blank card must appear. 

User-Prepared Subprograms 

In addition to preparing the input data cards described in the previous section, the 
user must prepare a group of subroutines and functions which specify the geometry and 
the load distribution. The structure of each individual subprogram is described sub- 
sequently in detail. 

In order to allow the use of parameters in the subprograms which may change from 
problem to problem in a sequence of problems a block of data called CONST (100), an 
array containing 100 elements, has been set aside. This array can be referenced in any 
subprogram by a COMMON statement. The locations of this array and others in 
COMMON which are required in preparing the subprograms are given in table 3. A 
description of the subprograms follows. Examples are presented in a later section. 

Geometry of the middle surface SUBROUTINE INPUT (NMAX).- The purpose of 
this subprogram is to calculate values for the following variables at the edges O and L 
and at each interior shell station; R, GAM, OMT, OMXI, and DEOMX where 


R(K) 


GAM(K) 


OMT(K) 

~ Mk 

OMXE(K) 

~ Mk 

DEOMX(K) 

" (“V)k 
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and where K is the shell station. In addition, the arc length increment DEL must be 
calculated as follows 

DEL " a(NMAX-2) 

where S is the maximum arc length. This information is transmitted to the main 
program by a COMMON statement. The location of these variables in COMMON is 
given in table 3. 

Any data required to calculate the preceding geometric quantities can be read by 
this subprogram or be transmitted from the main program through the CONST array. 

The quantity NMAX, the number of shell stations, is referenced through the argument 
list. 

Shell wall thickness and ratio of wall thickness to cover-plate thickness .- The shell 
wall considered in this program is as shown in figure 5 and discussed in appendix C. 

It is composed of three layers and is symmetric about the middle surface. It is assumed 
that the core is rigid normal to the middle surface but that it offers no resistance to 
extension or bending tangent to the middle surface. Four FUNCTION subprograms are 
to be specified which describe the shell wall construction so that the bending and mem- 
brane stiffnesses and the thermal loads may be calculated. The quantities to be specified 
by FUNCTION subprograms at station K are as follows: 


^ - - FUNCTION HHT(K,DEL) 

ho 

/— Y - FUNCTION DHHT(K,DEL) 

\ h o ) 

£ - FUNCTION HRA(K,DEL) 



- FUNCTION DHRA(K,DEL) 


where K is an index indicating the shell station and DEL is the arc length instrument A. 
For an isotropic shell of uniform thickness these quantities have the following values: 
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Variable parametric data can be transmitted to these functions by reading data into 
COMMON location BL14 (CONST(IOO)) in the main program and then transmitting the 
data to the subprograms by using a suitable COMMON statement. 

Temperature distribution.- The temperature, for any Fourier index, is assumed to 
be given in the following form: 

T (n) _ T (n) + AT (n)^ ( 32 ) 

where T^ is the temperature of the reference surface and AT^ is the temperature 
difference between inner and outer surfaces per unit thickness. This specification for 
this temperature distribution results in thermal forces and moments as presented in 
appendix D. The quantities T^ and AT^, and their derivatives with respect to £, 
are to be specified at the station K by functions as follows: 

tJ“) - FUNCTION TEMP(K,DEL) 

AT} 11 ) - FUNCTION DELT(K,DEL) 

T^ - FUNCTION DTEMP(K,DEL) 

AT^ - FUNCTION DDELT(K,DEL) 

where K and DEL are the same as defined previously. 

Surface loads.- The surface loads p^, p^, and p^ are specified at the sta- 
tion K by FUNCTION subprograms as follows: 

p - FUNCTION P(K,DEL) 
p^ - FUNCTION PX(K,DEL) 

P e - FUNCTION PT(K,DEL) 

where K and DEL are as defined in the preceding section. 

Program Output 

The output for this program is divided into three sections. The first section is 
devoted to a statement of the problem being run. It contains much of the input informa- 
tion and serves as a check on the accuracy of such information. The entire problem is 
summarized and the user may easily check the input data for errors. The data pre- 
sented in this section should be scrutinized before proceeding to an evaluation of the 
program calculations. 
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The second section presents a listing of shell geometrical data and is printed at 
the user's option. This output option is controlled by IND1 where geometrical data are 
printed if ESTDl = 0. The third section presents the data calculated by the program in 
nondimensional form. The relation between the nondimensional quantities and physical 
quantities for the output variables is given in table 4. Also given in table 4 is the title 
of the column heading under which the data appear on the output list and the location of 
each of the variables in storage at output time. The form of the output is presented in 
table 5. 


Errors 

The results should be checked to make sure that necessary conditions are satisfied. 
For instance, a simple check can usually be made to determine that an arbitrary portion 
of the shell is in a state of overall equilibrium. If the answers appear unreasonable, 
further checks should be made. 

Inconsistency of input data .- It is possible that the boundary conditions are incor- 
rectly specified or that incorrect input functions specifying geometry, loads, or tempera- 
ture are in use. Careful checking of output data describing shell parameters will help to 
eliminate many problems of this type. 

Numerical error.- There are two types of numerical error which are of importance. 
The first is the numerical error which is inherent in the size of the increment. The 
meridional increment must be sufficiently small so that a number of shell stations exist 
in any boundary layer associated with boundary condition effects. Another type of error, 
inherent in digital calculations, is the loss of significant digits due to truncation error as 
the number of shell stations increases. 

Estimation of Increment Size 

The practical problem of round-off error which occurs in the floating point arith- 
metic operations on a digital computer requires that attention be given to the estimation 
of a minimum increment size. Because of the complexity of the general problem it is 
not possible to write down an equation which relates minimum increment size to the 
geometric properties of the shell. 

A primary source of round-off error can occur in the calculation of the A, B, 
and C matrices, as defined by equations (20). It is seen that the elements of B, for 
instance, are formed by addition of appropriate elements of the E and G matrices 
(see eq. (13)). Furthermore, since the first term of the sum involves the element of E 
times the large number and the second term involves the element of G times the 
small number A, it is conceivable that for very small A the second term would have 
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no significance in the sum. Thus, for a homogeneous conical shell of constant thickness 
under axisymmetric loading, for example, it can be shown that the first element of the 
first row of B is of the following order: 

0(M) = + 0(A) (33) 

1 /l \ 

and since A « 1 the quantity — is of order Of — 1. In order for the second term on 

the right-hand side of equation (33) to have j significant figures in a sum where the 
first term has m significant figures, A can be estimated as follows: 

A 2 s* 10^' m 

Where m = 8 and j is desired to be 3, A is found to be of the order of 

-3 

A = 3 x 10 . Thus, for a very short cone, that is, a cone having a maximum nondimen - 

sional meridional length which is only one order of magnitude greater than the order of 
the minimum value of A calculated, only a relatively few stations can be considered. 

In order to check this rationale a very short truncated cone was considered, having 
a nose half-angle of £ radians and a nondimensional length approximately 20 times the 
calculated minimum increment length. The nondimensional radius to one edge was taken 
to be 1, whereas the radius to the other edge was specified to be 0.94545. This gave a 
nondimensional arc length of 0.0628. The edge having the smaller radius was free of 
forces; at the other edge the following conditions were prescribed: 

W R =1 

w A = 0 

= 0 

u 0 = o 

where Wr and are the radial and axial components of displacement, respectively. 

Since no axial forces are applied to the system, the axial component of the edge forces at 
the edge p = 1 should be zero. The axial-force resultant calculated for each increment 
is plotted against the increment size in figure 8. 

It can be seen that the residual axial force increases rapidly with decrease in 
increment size. For this particular case the largest increment which is the 
same order as the minimum increment length calculated in this example gives the small- 
est axial-force resultant. 

The message here is that, for some problems, increasing the number of stations 
does not increase the accuracy of the solution. In fact, for specific cases, a relatively 
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large error may result from using a large number of stations and thereby a small incre- 
ment size. 


General Considerations 


Suggested procedure .- It is suggested that when the INPUT subroutine is fairly com- 
plicated it be debugged separately until the user is convinced that the geometry is being 
correctly calculated. The program can then be easily converted to a subroutine. Simple 
problems should be considered first in order to gain familiarity with the potentialities of 
the program. 


P rogram running times .- On the average, a problem having 502 stations and having 
nonsimple geometry takes approximately 30 seconds on an IBM 7040/7094 II Direct 
Coupled System for one Fourier component. 


Program limitati ons and restrictions .- Most of these limitations have been men- 
tioned in various places throughout the text; however, they are summarized here for 
convenience. 


a. Young's modulus, Poisson's ratio, and the thermal expansion parameter 
are constant with respect to tangents to the middle surface. 



b. Three layers, symmetric with respect to the middle surface, are allowed. 

c. Only load and temperature distributions which can be represented as a Fourier 
series with respect to the circumferential coordinate can be considered. 


EXAMPLE PROBLEMS 


The example problems are typical of those which the user of the program might 
encounter in actual practice. The problems are meant to illustrate the preparation of 
input data, input subroutines, and boundary conditions. 


Cylindrical Shell With Edge Moments 

The problem considered here is that of a cylindrical shell of uniform thickness h Q , 


-2 


Eh, 


subjected to the end moments M fc = 10 

* 1 - 1,2 


cos n0 with Ut = Ug = W = 0 at each 


boundary. The calculations were made for — =50, v = 0.3, — = 1, where a is taken 

h 0 a 

to be the radius of the cylinder and S is the total arc length; 500 increments were used. 

Nondimensional parameters .- It is seen from table 4 that the relation between the 
physical and the nondimensional moment is given by 
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a o h o 3 


m, 


cos n0 


The value of is taken to be unity at the boundary for each Fourier index to give the 

following expression for ct 0 : 


10~ 2 aE 

h 0 (l - ^ 2 ) 


In general, the reference stress a Q is determined by an examination of the relationship 
between the dimensional and the dimensionless forcing function for a specific problem. 
The forcing function may be a surface pressure, a prescribed nonzero edge force or dis- 
placement, or a prescribed temperature variation. 

Boundary conditions.- For each value of the Fourier index the boundary conditions 
at !j= 0, 1 are m^ = 1, u^ = u^ = w = 0. The boundary condition matrices 12, A, 
and l are identical at each end and are as follows: 


n = [o] 


a = [i] 

l = [ 0 , 0 , 0 , 1 ] 

where I is the unit matrix. 


User -prepared subprogram s . - 

a. INPUT: At each station the following relations hold: 
p' = = 0. Also, the arc length increment is given by 

INPUT subroutine is as follows: 



co e = a; 

S 

a(N - 2)' 


A suitable 


SUBROUTINE INPUT (NMAX) 

COMMON R(502)/BL1/GAM(502),OMT(502),OMXI(502),DEOMX(502) 
1/BL3/NU, LAM, N,EALSIG, CHAR, DEL 
READ(5,1)SMAX 
1 FORMAT (F10. 5) 

DEL=SMAX/ FLOAT (NMAX - 2) 

DO 2 1=1, NMAX 
R(I) = 1. 

OMT(I)=l. 
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I 


OMXI(I)=0. 

DEOMX(I)=0. 

2 GAM(I)=0. 

RETURN 

END 

b. Thickness, thickness variation, ratio of total thickness to cover -plate thickness, 
and variation of ratio of total thickness to cover-plate thickness: A suitable set of func- 
tions to define these variables for a shell of uniform thickness h 0 is as follows: 

FUNCTION HHT(K,DEL) 

HHT=1. 

RETURN 

END 

FUNCTION DHHT(K,DEL) 

DHHT=0. 

RETURN 

END 

FUNCTION HRA(K,DEL) 

HRA=2. 

RETURN 

END 

FUNCTION DHRA(K,DEL) 

DHRA=0. 

RETURN 

END 

c. Surface loads: For the problem considered, p^ = p^ = p^ = 0. A suitable 

set of functions specifying these variables is as follows: 

FUNCTION P(K,DEL) 

P=0. 

RETURN 

END 
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FUNCTION PT(K,DEL) 

PT=0. 

RETURN 

END 

FUNCTION PX(K,DEL) 

PX=0. 

RETURN 

END 

d. Temperature: The temperature is uniform and is equal to a reference tempera- 
ture level at which no thermal strains are produced. A suitable set of functions which 
specify the temperature variables is as follows: 

FUNCTION TEMP(K,DEL) 

TEMP=0. 

RETURN 

END 

FUNCTION DTEMP(K,DEL) 

DTEMP=0. 

RETURN 

END 

FUNCTION DELT(K,DEL) 

DELT=0. 

RETURN 

END 

FUNCTION DDELT(K,DEL) 

DDELT=0. 

RETURN 

END 

Results.- The results of this calculation were compared with those presented in 
reference 1 for this shell configuration. The agreement was excellent. 
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Closed Spherical Shell Segment 

Membrane solution .- The problem shown in figure 9 is a closed spherical shell of 
uniform thickness h Q , subjected to a uniform hydrostatic external pressure q over the 
surface, and having the edge conditions = 0 which correspond to 

the membrane solution. The calculations were made for = 1000, v - 0.3, <p 0 = % 

ho 6 

and n = 0, where a is the shell radius; 50 increments were used. For this case, 
since the shell is closed at | = 0, the indicator IND5 is set equal to 1 and the finiteness 
conditions which are associated with the Fourier coefficient n = 0 are generated within 
the program. The results are in generally good agreement with the membrane solution 

03 , 

for the problem. The stress resultants and are each found to be equal to — 

in agreement with the membrane theory; however, a small residual moment on the order 

E O 

of 10 qh 0 is calculated, whereas the membrane theory is based on the assumption that 
the moment is identically equal to zero. 

Ela stic constraint at the boundary .- In the previous problem the edge conditions 
corresponded to those associated with the membrane solution. The present calculation 
is identical except for the boundary conditions. For this case the boundary conditions 
are taken in a more complex form in order to illustrate the specification of nondiagonal 
boundary condition matrices. The boundary conditions are taken to be 

Wy= 0 

= 0 

f H + kw H = 0 

u e = o 

where the components Wy, W H , and F H are defined in figure 9, and k is a propor- 
tionality constant. 

The horizontal and vertical components are related to shell variables as follows: 

F h = sin <p Q + cos <£ 0 
W H = W sin cp 0 + U| cos 4> 0 
Wy = W cos 0 O - sin 0 O 

The boundary conditions in terms of the nondimensional shell variables are as 
follows: 

-u^ sin <p Q + w cos d> 0 - 0 
u 0 = ° 
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where 


? I sin <p Q + cos 0 O + k(u^ cos 0 O + w sin 0 O ^ = 0 
m ^ = 0 


k 


a k 
ho E o 


The boundary condition matrices for this case are as follows: 


n L = 


0 

0 

COS 0 q 


0 


0 

0 

0 

0 


0 

0 

sin 0 O 
0 


0 

0 

0 

0 


-sin 0 O 0 

0 1 

Ajj = - 

k cos 0 O 0 

0 0 


COS 00 
0 

k sin 0 O 
0 


0 

0 

0 

1 





0 

0, 


The results N q and for the calculation of the sphere with hydrostatic pressure for 

these boundary conditions for k = 0.1 are presented in figures 10 and 11, respectively. 
These results show that the membrane results are reasonably valid through ^ = 0.8; 

O 

beyond that point a large boundary-layer effect is present, especially for N# and M^. 
An output listing for this problem for 50 equal increments is presented in appendix E. 


Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., January 5, 1967, 
124-08-06-11-23. 
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APPENDIX A 


DETAILS OF ANALYTIC FORMULATION 


The detailed derivation of the set of equilibrium equations (13) is presented in ref- 
erence 1, and is repeated here to give a complete presentation of the basis for the com- 
puter program. In proceeding from equations (8), (9), and (10) it is necessary to define 
strain displacement relations, stress-strain relations, and to separate variables. The 
derivation of equations (13) is given in the following sections. 


Strain Displacement Relations 

The displacements and rotations are shown in the positive sense in figure 3. The 


relation between rotations and displacements are as follows: 


= U- — + 

« a B? £ £ 


$ 


If 1 9W . ,, TT \ 
- a^ pW + 0U0 j 




The membrane strains are given by the following equations: 

i/ 9u i 


£ ri\Tf + "{ w 






' £0 2a\P d0 

The bending distortions are given by 


_ 1 

-jf 




- 5(5 


K Z6 2a 


d ®6 T 

+ 


(Al) 


(A2) 




(A3) 


1 9$t 9$ fl 1 , x /i 8Ut 9U„ \ 

P^/ + “aT " r *e + 2 i:K ' w e)\p ^ ' ~W ' yU 7 
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Stress -Strain Relations 

By using the usual Kirchhoff hypothesis the stress-strain relations are written as 
follows: 


„ ( a s “ va e) 

+ £ _ - 1 - + 




E 


oT 


(Oo ~ Wt) 

E U +0lr 

, +r K - LL±A a 


(A4) 


where ? is the normal distance from the middle surface and T the temperature 
change may vary with £ as well as with £ and 0. The forces and moments are 
assumed to be related to the stresses by the following integrals: 


N ri 

* 

^ d? 

M ri 

''1 



N e'j 


M e = i 

V e d? 

> 

(A 5) 


r "i* d ? 

M Z0 = 

I d5 _ 




where the integrations extend through the total shell thickness. 

Equations (A4) and (A5) are combined to give the following relations: 


n £ - . J 

EoT d? 

+ “ 


J Ed? 

J Ed? 


^ EaT d£ 

j’ E « 

jEd? 

(i + 


jEd? 

j 


(A6) 
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and 


M t - vM a 

u--* 


j* ?EoT d? 
J ? 2 E d? J ? 2 E d? 

j* ?EoT d? 


- vM t 

«e = — i + 


J ? 2 E d? J ? 2 E d? 


(1 + i/)M 


K te ~ 


M 


j* ? 2 E d? 


(A7) 


where it has been assumed that the middle surface is defined so that the following rela- 
tion is satisfied: 


J ?E d? = 0 


(A8) 


Fourier Series Expansion 


The complete set of field equations for the 17 unknown quantities N^, Ng, N^, M^, 

M e , M^ e , U^, U 0 , W, cp £, <p 0 , e^, e 0 , e^, k^, k q , and /c^ are given by 

equations (8) to (10), (Al), (A2), (A3), (A6), and (A7). It is assumed that the surface loads 
and temperature distribution can be expanded into a Fourier series in terms of the cir- 
cumferential coordinate. The loads and temperature distribution are given by the fol- 
lowing equations: 


q = Yj p^(£) cos ne 

n=0 

h °° 

q { = ^ P ? n)( « >COS nS > 

n=0 f 


(A9) 


q 


, oo 

_ CT o h o V 

0 a Zj 

n=l 



(£)sin n0 


OO 

T = ^ T^(|,?)cos n0 
n=0 

where <r 0 is a reference stress level and h 0 is a reference thickness. 


(A10) 
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The independent variables are thus expanded in the following nondimensional form: 


oo 

N £ = °o h o y 4 n)c ° s 


n0 


n=0 


°° 

N e = CT o h o X tf ] cos nO \ 


n=0 


OO 

= °o h o y f fe sin ne 


n=l 


(All) 


Mi 


_ °oV £ 


3 °° 




n=0 


m^cos n0 


M e = ^a 2 - Z m e n)cos nd / 

n=0 


-T _ a oh 0 ° V 
he-~T~ L 


M>. „ = > ml n isinn0 




n=l 


(A12) 


U, = ^y uWcos ne 
I E 0 L % 
n=0 


U, 


OO 

E 0 L, 
n=l 


g^sin n0 


W = — Y w^cos ne 

En Li 


n=0 


(A13) 


e> = — y ej^cos ne 
£ E o L, Z 

° n=0 


: fl = — y el n ^cos ne 

0 Eo L 8 




n=0 


j 


(A14) 


(Equations continued on next page) 
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e 


4 e E r 



sin n0 


$ 


$ 


= fo 
4 E 0 

y 

i—j 

<^^cos nd 


n=0 


_ a o 
6 E 0 


sin n0 


n=l 

J 


(A 14) 


(A 15) 


K* = Y kl^cos n 9 
4 aE n L 4 
° n=0 

K e * I 4” )cos " e 1 

n=0 


*40"^ I 

n=l 



sin n0 


(A16) 


By using equations (All), (A12), and (A9), the equilibrium equations (8), (9), and 
(10) can be decoupled for each Fourier index n. Omitting the superscript (n), the set 
of partial differential equilibrium equations (8) to (10) become the following set of coupled 
ordinary differential equations which govern the Fourier coefficients: 


*{’ + ' *«) + (l)‘«e + x2 [“{ m {' + - m e) + (fe)( 3 " { - "e) m se] 

V + 2 yt i$ - (?)*e + x 2 {-(|) a V n e + |( 3 "e - "jJ’V 

* |[’ , ( 3 "e + "{) ‘ “«] m {<} + Pe - 0 


+ pi = o'! 


> (A 17) 




0 J e t e 


! { r 


+ + 2Xm 






- (a^u) e m^ + 



‘ yn V + (ir) m |e' + (^p") m se} + p - 0 


j 
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where 


X = 


(A18) 


Similarly, using equations (A13), (A14), (A15), and (A16), equations (Al), (A2), 
and (A3) yield the following additional equations for each Fourier coefficient: 



cp Q = I -)w + w fl u. 


(A 19) 


e fc = u fc ' w fc w 

? ? ? 


e e = [p) u e + ^ + V* 


e i« " I 1 


V ■ e - [prs 


(A20) 


k « = 


k e = (ye + ^ 


k - - 

he~ 2 




-$*•( + *» - ™><, + 1 ("9 - “ £ ) (j “ 4 + V + ru| 

J 


The Fourier coefficients for forces and moments are found to be 


and 


t^ = - t{^ 


t e = b ( e e + ^i) 
Ue = b(1 " u)e ^e 


(n) 


(A21) 


(A22) 


j 


m^ = d^k^ + vk ^ - m^ 
m 0 = d^kg + uk^j - m^ 

m^ e = d(l - v)* ke 


c\ 


(A23) 
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j* ? 2 E d? 
E o h o 3 ( 1 " v 2 ) 


(A24) 


(A25) 


T ^oV 1 “ ") 


(A26) 


m (n) . 

T ' "A’d-i 


(A27) 


The set of 17 field equations for each Fourier number n for the Fourier coeffi- 
cients tg, t 0 , tgg, m^, m 0 , m^g, u^, u 0 , w, <p Q , e^, e 0 , e^g, k^, kg, 

k^g is given by the 17 equations (A17) and (A19) to (A23). 


Reduction to Set of Second-Order Differential Equations 

The set of 17 field equations constitutes an eighth-order system. By combining 
equations to eliminate variables, the 17 equations can be reduced to the following set of 
4 second-order differential equations in terms of the variables u^, Ug, w, and m^: 

Ez" + Fz’ + Gz = e (A28) 

where the relation 

mg = + d(l - ^ 2 )kg - (1 - (A29) 

has been used and where 


(A30) 
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and 



Eh 

0 

0 

0 


0 


e 22 

E 23 

0 


0 


e 32 

E 33 

E 34 


0 


0 

e 43 

0 


F 11 

F 12 

F 13 

F 14 


f 21 

F 22 

F 23 

0 


f 31 

F 32 

F 33 

F 34 


f 41 

0 

f 43 

0 


G 11 

g 12 

G 13 

g 14 


g 21 

g 22 

g 23 

g 24 


g 31 

g 32 

G 33 

G 34 


g 41 

G 42 

G 43 

g 44 


" e r 





= < 

e 2 

u 

> 





e 4 






(A31) 


The nonzero elements of the E, F, G, and e matrices are defined in appendix A of 
reference 1. 




Boundary Conditions 

In the Sander theory the expressions for virtual work per unit length at boundaries 
^max are 


W B = ± (N^ + N^ 0 U 0 + Q^W + 
where Wg is the virtual work per unit length and 


(A3 2) 
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*«s- N «e + [(^)-(5q3 M 


l te 


and 


Qi 


‘{k)[kY pM i) + 2 K^ ] - p ' Ml 


(A3 3) 


(A34) 


Expanding and in Fourier series and normalizing gives 


*40 = *40 + 


( t )( 3 “8 • "{)“{« 


- >2(~ 
= a n 


m^' + y + 




(A3 5) 


The forces defined by equations (A33) and (A34) are effective inplane force and 
transverse shears per unit length, respectively. The forces and moments acting on the 
boundary are shown in figure 4. This form of the virtual work of the forces indicates that 
either the effective force (moment) or related displacement (rotation) may be prescribed 
or, more generally, the force and displacements may be linearly related. The linear 
relation provides a means of handling elastic restraints. A convenient boundary condi- 
tion representation expressed in nondimensional form can thus be written for the nth 
Fourier coefficient as follows: 

f2y + Az = l (A36) 


where y is the force vector 


y 

a 


(A3 7) 


The vector y is related to z 


The nonzero coefficients of H, 


as follows 

y = Hz' + Jz + f (A38) 

f, and J are presented in appendix A of reference 1. 
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FORTRAN IV SOURCE DECK 


The listing of the source decks for the main program is on the following pages 



appendix b 


C MAIN PROGRAM 

C THIS PROGRAM CONSISTS UP THE MAIN PROGRAM TOGETHER WITH THE FOLLOWING 

C SUBROUTINES 

C MAT I NV 

C EFG 

C OUTPUT 

C STRESS 

C EQ7 3 

C INIT 

C FINAL 

C FORCE 

C PANOX 

C A8CG 

C HFJ 

C KLT 

C BOB 

C POLE 

C IN ADDITION TO THE ABOVE THt USER MUST SUPPLY THE FOLLOWING SUB 

C RGUTINtS AND FUNCTIUNS. 

C SUBROUTINE INPUT — CALCULATES THE SHELL GEOMETRY* 

C FUNCTION P— SPEC I FILS THE NORMAL PRESSURE DISTRIBUTION* 

C FUNCTION PX— SPECIFIES THt MERIDIONAL PRESSURE DISTRIBUTION. 

C FUNCTION PT-SPECIFIES THE CIRCUMFERENTIAL PRESSURE DISTRIBUTION 

C FUNCTION TEMP-SPEcl FI ES THE MERIDIONAL TEMPERATURE DISTRIBUTION 

C FUNCTION DTE MP— SPECIFIES IHE DERIVATIVE OF THE TEMPERATURE 

C FUNCTION DEL T- SPECIFIES THE DISTRIBUTION OF THF TEMPERATURE 

C VARIATION THROUGH THE THICKNESS. 

C FUNCTION DDFLT- SPECIFIES RATE OF CHANGE OF DELT 

C FUNCTION HHT-SP EC i F I ES THE SHF LL THICKNESS DISTRIBUTION 

C FUNCTION DHHT-SPECIFitS THE DERIVATIVE OF TFE THICKNESS. 

C FUNCTION HR A-S P EC I F I ES THE DISTRIBUTION OF THF RATIO OF THF 

C TOTAL THICKNESS- TO -COVER PLATE THICKNESS. 

C FUNCTION DHR A— THE DERIVATIVE OF THE ABOVE RATIO. 

INTEGER FREO 
REAL NU , L AM , N # J A Y 

COMMON R(502)/GL1/GAM(502) ,UMT(502) , OMXI (502) ,DE0MX(502) 

1/BL2/E(4 f 4) ,G(4,4) ,FU, 4) 

2/BL3/NU, LAM,N, FALSIG,CHAR f DEL 
3/BL4/CE F ( 4 ) 

4/8L5/H(4,4),EF(4),JAY(4,4) 

5/BLo/A(4 f 4),B(4,4),t{4,4) , SMAG(4) 

6/BL7/PFE{4,4 T 502)»X(4,bG2) 

7/BLV/OMEG 1( 4,4) ,CAPL1(4 ? 4) # EL1(4) 

8/8L11/0MEGL ( 4 , 4 ) , L APLL ( 4 , 4 ) , EL L < 4 ) 

9/BL13/AK( 3*4) , ALL (3,4) , STH ER( 3 ) /BL 14 /CON ST ( 100 ) 

DIMENSION Z(4,50?) , IPIVOT (4), I N0EXI4,2) 

EQUIVALENCE! X( I, 1) ) 

1 FORMAT ( 5F6.3»6I4) 

2 FURMATt 1H1, 16H PROBLEM NUMBER 16///) 

21 FORMAT (4E16.8) 

31 FORMAT ( 72H 

1 ) 

50 FORMAT (214) 

60 FORMAT! 1H0 , 5 5X 10H INPUT DATA///) 

100 FORMAT ( 52 H THE FOLLOWING CONSTANTS APPEAR IN COMMON BLOCK BL14///( 
11X1 P8F 16*8) ) 
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164 FORMAT ( IX 2CH REFERENCE LENGTH =1PE16. 7//1X 20H POISSONS RATIO 
11PE16.7// 1X20H TEMP COEFFICIENT =1PE16. 7// 1X20H REF. THICKNESS 

2 =1PE16.7//IX20H FOURIER INDEX = 1PE 16 . 7 // 1 X20H NUMBER OF ST ATI 

30NS=I 5/// ) 

165 FORMAT ( 34H THE INDICATORS ARE SET AS FOLLOW S//7H IND1=I4,3X7H IN 
102=14, 3X7H IN03=I4,3X7H IND4=I4,3X7H IN05=I4//) 

12 RE Al)( 5 , 1 ) NU»TKN,CHAR,N,EALSIG,N0,INDl,IND2»IND3»IND4,IN05 
WR I TE ( 6 » 2 ) NO 
READ (5,31) 

WR I T E ( 6 , 3 1 ) 

I F ( IN02.EQ.0) GO TO 70 

READ! 5,21 ) (CONST! I ) , 1=1, IN02) 

70 IF ( IND3.NF.O)GO TO 3 
GO TO 4 

3 REAO( 5, 50 )NMAX,FREW 
CALL INPUT(NMAX) 

4 LAM=TKN/CHAR 
WRITE (6,60) 

WRITE (6,1 64 )CHAR,NU, EALSIG, TKN , N , NMAX 

WRITE (6, 165) IN01, 1ND2, I ND3 , I N04, I ND5 
IF ( IN02.NF.0) 

1WRITE(6, 100) (CONST l i ),1=1,IND2) 

CALL HFJ( l, I NO 5, NMAX, 2.) 

CALL EFG( 2, I NO 5, NMAX ) 

CALL FORCE! 2, IN05.NMAX) 

CALL A3CG 

CALL I N I T ( IND5, IND4) 

NMAX1 =NMAX— 1 

DU 6 K = 3 , NM AX 1 

CALL EFG(K, I NO 5, NMAX ) 

CALL FORC E(K, IND5,NMAX) 

CALL A3CG 

6 CALL PANDX(K) 

CALL HFJ(NMAX,I ND5 , NMAX, 2.) 

CALL FINAUNMAX, INU5, INU4) 

DO 7 L = 2 , NMAX 1 
K=NMAX+ 1-L 

7 CALL E073 ( K ) 

IF! IND5.EQ.0.0R.IN05.EQ.3) GU TO 14 
CALL E Q73 ( 1 ) 

GO TO 15 

14 CALL EFG( 2,IND5,NMAX) 

CALL FORCE! 2, IND5, NMAX) 

CALL ABCG 

DO a 1=1,4 

51 = 3. 

52 = 0. 

DO 9 J= 1 , 4 

S1=S l+A( I , J)*Z( J,3) 

9 S2=S2+B(I ,J)*Z(J,2) 

8 SMAG( I )=SMAG( I )-Sl-S2 

CALL M ATI N V (C, 4, SMAG,1, DETERM, IPI VOT, INDEX , 4, I SCALE) 

DO 10 1=1,4 
10 Z( I , 1) =SMAG( I ) 

15 CALL STRESS! FREQ, NMAX, N, IND5) 

CALL OUTPUT ! FREQ , NMAX , I NUl , DEL , NO ) 

GO TO 12 

END 
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SUB ROUT INF HFJ(K, IND5.NMAX ,YAH ) 

C SUBROUTINE HFJ THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE H,F, 

C AND J MATRICES, AS DEFINED IN APPENDIX A OF REFE RENCE ( 1 ) , AT THE STATION 
C SPECIFIED BY THE INDEX K. 

COMMON R( 502) /BL1/GAMI502) ,GMT ( SC 2) , OMX II 50 2) ,DEOMX( 502) /BL3/NU, 
1LAM,N,EALSIG,CHAR,0EL/BL5/H14,4),FF (4) ,JAY{ 4,4) 

REAL NU*LAM,N, JAY, L2 

CALL BOB ( K , DEL , NU ,B,OB,D,DD) 

L2=LAM**2 
D 1= ( l.-NU) 

OX=OMX T ( K ) 

REG=0 . 

IF ( YAH . EQ. 2 . ) REG= 1 , 

EAL=EALSIG 
T = T EMP ( K* DEL ) 

DLT=DELT I K » DEL ) 

H1=HHT I K* DEL ) 

HRB=HR A I K » DEL ) 

FF ( l)=-2.*Hl*EAL/ (Dl*HRb)*T 
FF ( 2 ) =0 . 

FF(3)=L2*CHAR*EAL*uLT*Hl**3/3.*(1.5/HRB-3./FRB**2+2./HRBt‘*3) 

FF ( 4 ) =0 . 

I F ( IND5 ) 1 ,1,2 

2 IF ( ( ( IND5-2) .LF .0) . AND. CK.EU. 1 ) ) GO TO 8 
I F ( ( ( IN05-2) .GT.O.AND. (K.EU.NMAX) ) GO TO 8 
GA=GAM1K) 

FF ( 3 ) =FF ( 3 ) *GA 
I RA=R ( K ) 

OT=OMT ( K ) 

ENR=N/RA 

0XT = 3. #UMX I I K ) — CMT ( K ) 

0TX=3.*0MT(K)— OMXIIK.) 

0L=D*L2*D1*FNR 
H( I , I) =B 
H( 1 , 2 ) =0. 

HI 1 , 3 ) =0. 

HI l,4)=0. 

H I 2 , I ) =0 . 

H(2,2)=B*DI/2.+L2< t O*Dl/8.*UTX**2*REG 
HI 2, 3>=DL/?.*0TX*RcG 
HI 2 , 4 ) =0. 

HI3, 1)=0. 

HI 3,2) = OL *OTX*Y AH/4 . 

FNR2=ENR* *2 

HI 3,3)=L2*D*DI*( YAH*ENR2+( I . + NU ) *GA**2 ) 

GA2=GA**2 
H(3,4)=L2 
H I 4 , 1 ) = 0 . 

H(4, 2 ) =0. 

H(4,3)=-l. 

H(4 , 4 ) =0. 

JAY! 1, 1)=NU*GA*B 
JAY I 1,2) =NU*B*ENR 
JAYI l,3) = B*(OX«-NU*UT) 

JAY! 1 , 4 ) =0. 

JAY l 2, 1 )=-B*Dl*ENK/2.-DL/8.*UXT*OTX*REG 
JAY I 2 , 2 ) =-GA*H( 2,2) 

JAY ( 2 , 3 ) =-GA*H( 2,3) 

JAY ( 2, 4 ) =0 . 

JAY I 3, 1 ) =— L 2*D*D 1* I I 1 . + NU) *GA2*OX+ENR2/4. *0 XT* YAH) 

JAY! 3»2) = -GA*DL/2.*( 2.*OT* < 1.+NUJ+OTX/2 .*Y AH ) 

JAY I 3,3)=— L2*D*DI*( 1 . +NU + YAH ) *GA*ENR2 
JAY! 3,4) =L 2*D I*GA 


JAYJ4, 1)=0X 
JAY ( 4 » 2 ) = G • 

JAY(4,3)=0. 

JAY(4»4)=0. 

GO TO 5 
8 DO 61=1,4 
DO 6 J=1 , 4 
H ( I,J)=0. 

6 JAY(I,J)=0. 

H ( 1, 1 ) = B* ( 1 . + NU ) 

H ( 1,2) =N* 8*NU 

H(2, l ) = B*0 1* ( -N)/2. 

AW3=Dl*<-3.*{ l.+NU)+N**2*(3.+NU) ) 

AW=AWB*D 

Cl=AW/(-2.+NU*( N **2—1 . ) ) 

ATHETA=1. 5*0*N*0X*D1*(3.+NU) 

AX 1= 1 • 5*D*OX*Dl* ( 2 • * ( 1 . +NU ) -N ) 

H(3»4)=L2*( 2.— NU+2 • *AwB ) 

H(4, 3) =-l . 

J AY ( 1 » 3 ) = B* ( 1 . + NU ) *UX 
JAY( 4, 1 )=OX 
OH=OHHT(M,OEL) 

DTOH=OHR A ( M , DEL > 

DDLT=DDELT(M,DEL) 

DTMGN=CHAR*EALSIG/3./0I*( ( 1.5/HRB-3./HRB**2+2./HRB**3)*( H1**3*0DL T 
1 + 3.*DH*H1**2*DLT)+l)L T*H1**3*0T0H/HRB**2*(-1 .5+6. /HRB— 6. / 

2HR8**2) ) 

FF ( 3 ) = 0 T M DM * ( 2 • * A w B / (-2 . + NU* ( N**2-l . ) ) +2.-NU) 

5 RETURN 
END 
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SUBROUTINE: EFG( K, 1N05,NMAX) 

C SUBROUTINE EFG THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE E * F , AND G 
C MATRICES, AS DEFINED IN APPENDIX A OF REFERENCE (1),AT THE STATION SPECIFIED 
C BY THE INDEX K. 

COMMON R ( 502) /BL1/GAM( 502) ,UMT (502) ,OMXI ( 50 2) ,DEOMX( 502) /BL2/E(4,4 
I) , G ( 4 , 4 ) ,F(4»4) /BL3/NU,LAM,N,EALSIG,CHAR,0EL 
REAL NU , LAM, N , L AM2 
CALL BDB( K,DFL,NU,B,OB,D,DD) 

E(l»i)=B 
E ( I , 2) =0. 

E ( I , 3 ) =0. 

E ( 1 , 4 ) =0 . 

E ( 2 , I ) =0. 

Dl=( l.-NU) 

LAM2=L AM**2 
RA=R ( K ) 

GA=GAM(K) 

OX = OMX I ( K. ) 

OT=OMT ( K ) 

DEX=DEQMX ( K ) 

GA2=GA**2 
REX=( 3.*OT-OX) 

RXE= ( 3 .#OX— OT ) 

OTX = OT *OX 

DNLR=L AM2*D*N*Dl/( 2.*RA) 

DDNLR = l)NLR*OD/D 

E ( 2,2)=B*Dl/2.+LAM2*U*Dl*R£X**2/8. 

E(2,3)=QNLR*PEX 

£(2,4)=0. 

E ( 3 , l ) =0. 

E ( 3 » 2 ) =F ( 2,3) 

RAN= ( N/RA ) 

E ( 3, 3) =LAM2*0#D1*( 2. *KAN+ (l.+NU)* GA 2) 

E(3,4)=LAM? 

E ( 4 , 1)=0. 

E ( 4 , 2 ) =0. 

E ( 4 , 3 ) =-0 
F(4,4)=0. 

F ( I , I ) =GA*B+DB 

F ( I , 2 ) = ( 1 . +NU ) *R*N/ ( 2.*RA ) + DNL P*REX*RXE/4. 

F ( 1 , 3 ) =8* (OX+NU*OT )+LAM2*D*Ul* ( ( 1 . + NU) *GA2*CX+RAN*RXE/2. ) 

F ( 1 , 4 ) =L AM2*0X 
F ( 2 , I ) =— F ( 1,2) 

F(2,2) = (f)l/2. )*(GA*B+OB)-(LAM2*D*Dl*RFX/8.) *( 2 . *DE X-G A* ( 5 . *0 X 
1-3. *UT) ) + LAM?*0!)*01*REX**2/ci. 

F ( 2 , 3) =DNLR* ( 2. *( l.+NU) *GA*uT-DEX + 3 . *GA*( OX-OT > ) +DDNLR*R EX 
F ( 2 , 4 ) =0 . 

F ( 3 , 1 > =— F (1,3) 

F( 3, 2 ) =ONLR*( 3. *GA*OX-oA*uT* ( 5 .+2 . *NU) -DEX ) +DDNL R*REX 
F( 3, 3>=-LAM2*D*Dl*< ( 1 . + NU ) * ( 2 . *GA *OX*OT+GA* *3 ) +2.*GA*RAN) 

1+LAM2*DD*D1*( ( 1 ,+NU ) *GA 2+2 . *R AN ) 

F( 3,4)=LAM2*GA*( 2.-NU) 

F ( 4 , 1>=0*0X 
F(4,2)=0. 

F(4, 3)=— D*NU*GA 
F ( 4 , 4 ) =0. 

G( 1, 1 )=NU*0B*GA-NU*B*0TX-B*GA2-01*8*RAN/2.-LAM2*D*Dl*{ ( l.+NU )*GA2* 
10X**2+RXF**2*RAN/6. ) 

G( 1 , 2 ) =NU*N*Oft/RA-( 3.-NU ) / ( 2 . *RA) *GA*B*N-DN LR*2 . *GA* ( REX*RXE/ 8. 
l+( l.+NU ) *OTX ) 

G( 1 , 3 ) =B* { DEX+GA*( UX— UT ) )+D6*( GX+NU*OT )-L AM 2*D*D1*GA*RAN* (RXE/2.+ ( 
ll.+NU)*OX) 

G( l,4)=LAM2*Di*GA*0X 
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G( 2, l)=-3*GA*N*< 3.-NU)/( 2 . *RA > -D1*N*DB /< 2 . * PA ) +DNLR*2 . * ( - 1.* ( 1.+ 
1NU) *GA*0TX+GA/3.*< 6. *UTX-7 . *UX**2- 3. *0T**2 > -DEX/4.*( 5.*0T-3. *0X) ) 
2-DDNLR/4.*REX : *RXE 

G(2, 2>=-GA*F(2, 2 ) +U 1/2 . *a*0TX-B*RAN-LAM2*D* Cl* ( ( l.+NU)*0T**2*RAN 
l-GTX/8 . *R FX**2 ) 

G(2, 3)=-B*N*(OT+NU*UX)/RA+UNLR*(GA*OEX-2.*GA2*OX-2.*( i.+NU)*OT 
1 *RA,N + RE X* l GA2+0TX ) ) — L)UNLK*RfcX*GA 
G { 2 t 4 ) = -NU*L AM2*0T*N/RA 

G ( 3» l)=-B*GA*(0T+iNU*UX)+LAM2*0*Dl*JGA*( l. + NL>*(-GA*DEX+GA2*0X 

1- UX*RAN+2 . *OTX*OX ) +RAN/2 .*( GA*OX-GA*QT -3 .*DEX) ) 

2- LAM2*D0*Dl*( ( 1 . + NUJ*i J A2*OX+RAN/2.*RXE) 
G(3»2)=-B*N*(GT+NU*GX)/RA+UNLR*(2.*(l.+NU)*(OTX*OT-GA2*QX+2.*GA2 

1*0T-UT*RAN) +GA*DEX+3.*GA2*(UT-0X> +0T X* REX > - DDNLR* < 2.*< l.+NU) *GA 
2*0T +GA*RE X ) 

G( 3,3 >=-B*(OX**2+2.*NU*GTX+OT**2)+LAM2*D*01 *RAN*( ( l.+NU) *(OTX-RAN 
l+2.*GA2)+2.*(GA2+uTX ) ) -LAM2*DD*D 1 *R AN* ( 3 . + NU) *GA 
G(3,4)=-LAM2*(D1*GTX+NU*RAN) 

G ( 4 » 1)=D*(DEX+NU*GA*0X) 

G( 4, 2)=D*NU*N*0T/RA 
G(4» 3) =ll*NU*RAN 
G ( 4 . 4 ) =— 1 • 

RETURN 

END 


44 


non 


SUBROUTINE FORCE ( K, I ND5 , NMAX ) 

C SUBROUTINE FORCE THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE LOWER CASE 
C E-VECTOR AS DEFINED IN APPENDIX A OF REFERENCE! 1) , AT THE STATION SPECIFIED 
C BY THE INDEX K. 

COMMON R( 502) /BL 1/GAMl 502 ) ,OMT(502> ,OMXI (502) ,DEOMX( 502) /BL3/NU, 
1LAM,N,EALSIG,CHAR,0EL/BL4/CEE(4) 

REAL NUtLAMt N»L2 
RA=R ( K ) 

GA=GAM( K ) 

OX=OMX I ( K ) 

OT=OMT ( K ) 

T=TEMP( K» DEL ) 

DT=DTEMP( K,OEL) 

DELT1=DELT(K,DEL) 

DLT I=ODEL T( K,DEL ) 

PX1 = PX ( «.* DEL ) 

PT1=PT ( K* DEL ) 

P1=P(K,DEL) 

H=HHT ( K » DEL ) 

DH=DHHT ( K » DEL ) 

HRB = HR A ( K » DEL ) 

DHRB=DHRA ( K » DEL ) 

01= I NU 
L2=LAM**2 
EAL=EALSI G 

CEE(4)=CHAR*EAL*DELT1*H**3* 

1 (3./( 2.*HRB)-3./HRB**2+2./HRB**3) /3./D1 

TSUBT=2.*H*EAL/( D1*HRB»*T 

CEE ( l)=-PXl+2.*EAL*(H*DT+T*DH) /( D 1*HRB ) -TSU BT*DHRB/HRB 
CEE(2)=-PT1-(N/RA)*TSUBT-L2*D1*(N/RA)*0T*CEE<4) 

CEE( 3)=-Pl-(0X+0T)*TSUBT— L2*D1 *CH AR*EALS I G* ( ( 1 . 5/HRB-3. /HRB**2 
l+2./HRB**3) *(GA*H**3*ULTl+( 3.*H**2*GA*DH+H**3*(0X*0T-(N/RA)**2) ) ) 
2*DELTH-DELTl*H**3*DHRB 

3 /HRB**2*(-1.5+6./HRB-6./HRB**2> )/3./Dl 

RETURN 
END 


SUBROUTINEABCG 

SUBROUTINE ABCG — THIS SUBROUTINE CALCULATES THE A,B,C, AND LOWER 
CASE G MATRICES USING THE CURRENT VALUES OF THE E,F,G, AND LOWER CASE 
E MATRICES. 

COMMON /BL 4 /C EE (4) /BL2/E(4»4) »G (4,4) , F ( 4 , 4 > / BL6/ A ( 4,4 > ,B(4,4) , 

1C (4,4) , SMAG( 4) /BL 3/NU, LAM, N,EALSIG, CHAR, DEL 
REAL N » NU »L AM 
02=2. /DEL 
D4=4./DEL 
DX=2.*0EL 
0011=1,4 

SMAG(I )=DX*CFE( I ) 

D01 J= 1,4 

B ( I , J ) =-D4*E ( I , J ) +DX*G( I , J ) 

C( I , J ) =D2*E ( I , J )— F ( i,J) 

1 A(I,J)=D2*E( I » J ) +F ( I , J ) 

RETURN 

END 
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SUBROUTINE I N I T ( I NO 5 , IND4 ) 

CALCULATION OF THE P-MATRIX AND THE X-VECTOR AT THE FIRST STATION. 
COMMON/ BL 5/H(4»4I »FF (4) , JAY l 4, 41 /BL 7/PEE ( 4 , 4, 502 ) , X I 4 , 502 ) /8L9/ 
10MEG1(4,4) ,CAPL1(4,4) ,ELU,4> 

2/BL3/NU,LAM,N,EALSIG, CHAR, DEL 
3/BL2/E(4,4) ,G(4,4) ,F(4,4) 

4/BL4/CEEI 4) 

5/8L6/AI 4,4) , B ( 4, 4 ) , C ( 4 , 4 ) , SMAG ( 4) 

DIMENSION IPIV0TI4) , I NUEX (4, 2 ) , AO ! 4, 4) , 80 ( 4 ,4) , A3! 4, 4) , A4(4,4> ,G0( 
14) 

REAL J AY » N » NU » L AM 
21 F0RMATI4E 16.8) 

166 FORMAT! 1H0,46H THE BOUNDARY CONDITIONS AT S=0 ARE AS FOLLOWS//) 

167 FORMAT ( 21X5H0MEGA49X5HLAMDA34X3HELL///) 

169 FORMAT! IP 4E 1 2 . 4 , 6X 1P4E 12. 4, 6X 1 PE 1 2. 4 ) 

DO 8 1=1,4 
GO! I ) = 0 . 

EL1U)=0. 

DO 8 J = 1 , 4 
CAPL1! I , J ) =0 . 

8 0MEG1! I , J ) =0 . 

IF ! IND5 ) 9,9,10 

10 IF ( I ND5-2 ) 1 1 , 1 1 , 9 

11 CALL POLE ( N, DEL , F ,G) 

WR I TE ( 6 , 1 66 ) 

WR I TE ( 6 , 40 ) 

40 FORMAT! 1H0, 52H THE CONDITIONS FOR A SHELL POLE HAVE 8FEN GENERATE 
ID/// ) 

GO TO 16 

9 IFIIND4) 13,13,14 

13 R EAD I 5 , 2 I ) OMEGH 1 , 1 > , UMEG1 ! 2 , 2 ) , 0MEG1 ( 3 , 3 ) ,0MEG1< 4, 4 ) 

READ (5,21) C APL 1(1, 1 ) ,C APL 1 ( 2 , 2) , CARL 1 ( 3 , 3 > ,CAPL1!4,4) 

READ (5,21 ) ELI 

GO TO 15 

14 READ!5,2l )0MEG1»CAPL1»EL1 

15 WR I TE ( 6 , 166 ) 

WRITE(6, 167) 

DO 168 1=1,4 

168 WRITE16, 169)0MEG1! I , 1) ,0MEG1( 1,2) , OMEG 1 ( I , 3 ) ,OMEG 1 ( I , 4) , 

1 C APL 1 ( I, 1) , C APL 1 ( I ,2) ,CAPL 1(1,3 ) ,CAPL 1 (1,4) , 

2 EL 1 ( I ) 

DO 1 1=1,4 

DO 1 J = 1 , 4 

80! l,J)=H!I,J)/DEL+JAY(I,J)/2. 

1 AO! I,J)=JAY( I,J)/2.-H(I,J)/DfcL 
DO 2 1=1,4 

DO 2 J= 1,4 
S1 = 0. 

S2=0. 

DO 3 L= 1,4 

S1 = S1+0MEG1< I ,L)*A0(L, J ) 

3 S2=S2+0MEG1( I ,L)*BO(L,J) 

A3! I , J ) =Sl 

2 A4( I , J ) =S2 
DO 5 1=1,4 
Si = 0. 

DO 4 J= 1 , 4 

S1=S1+0MEG1( I »J)*FF(J) 

BO! I , J)=A3( I , J)+CAPL1( I, J)/2. 

4 AO! I , J)=A4( I, J) +CAPL1! I , J)/2. 

5 GO! I ) = EL1 ! I )-Sl 

CALL MAT I NV ( C , 4, CEE , 0,DE TERM, IPI VOT , INDEX , 4, I SCALE ) 

DO 50 1=1,4 



DO 50 J= 1,4 
S1=0. 

S2 = 0 • 

S3=0. 

DO 51 K=l,4 
S1=S1+C(I , K ) *R ( K , J ) 

S2=S2+C(I ,K)*A(K,J> 

51 S3=S3+B0( 1«K)*C<K« J) 

A3 ( I , J ) =S 1 

A4( I,J)=S2 
50 E ( I , J ) = S3 
DO 52 1=1,4 
DO 52 J = 1 » 4 

51 = 0. 

52 = 0. 

DO 53 K = 1 , 4 
S 1 = S 1+ B0[ I ,K)*A3[K,J) 

53 S2=S2+B0( I,K)*A4<K,J) 

F ( I , J ) =S1— AO C I , J) 

52 G( I , J ) = S2 

16 CALL MATINV(F,4,G,4fDfcTERM,lPI VOT , INDEX, 4, I SCALE) 
IFl IND5.FQ.0.0R.IND5.EU.3) GO TO 59 
DO 60 1=1,4 
XU , 1)=0. 

DO 60 J = 1 , 4 

60 PEE ( I , J, 1 >=G< I , J) 

CALL P AND X ( ? ) 

RETURN 

59 DO 54 1=1,4 
S1 = 0. 

DO 61 J = 1 , 4 

PEE ( I, J,2)=G< I , J) 

61 S1=E( I,J)*SMAG( J)+S1 

54 G0( I )=S1-G0< I ) 

DO 56 1=1,4 

S 1=0 • 

DO 57 J = 1 , 4 
57 S1=S1+F( I,J)*GO( J) 

56 X(I , 2)=S1 
RETURN 
END 
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SUBROUTINE PANOX(K) 

C SUBROUTINE PANOX THIS SUBROUTINE CALCULATES THE P MATRIX AND THE X-VECTOR 
C AT THE STATION SPECIFIED BY THE INOEX K USING EQUATION 129) OF THE TEXT. 
C0MM0N/BL6/A(4,4> ,B(4,4> ,C (4, 4 ) , SMAGI 4 ) /BL7 /PEE ( 4,4, 502) , X (4 ,502 ) 

01 ME NS I ON PI (4,4) » IPI VQT (4) , INDEX! 4, 2) , XI (4 ) ,X2I4) 

DO 1 1=1*4 
DO I J= 1 ? 4 
SUM=0. 

DO 2 L= l * 4 

2 SUM=SUM+C C I»L)*PEE(L» J »K— I ) 

I PI(I,J)=B(I*J )— SUM 

CALLMATINVI P 1 ,4 *X2* 0 ,DETERM, I P IVOT* INDEX,4» ISCALE) 

DO 4 1=1,4 
SUM=0. 

DO 3 J= 1 , 4 

3 SUM=SUM+Cf I ,J)*XC J,K-1) 

4 Xi(I)=SMAG(I )— SUM 
DO 5 1=1,4 

DO 5 J = 1 » 4 
SUM=0. 

DO 6 L= 1,4 

6 SUM=SUM+P 1 ( I , L ) *A ( L , J ) 

5 PEE ( I , J ,K )=SUM 
DO 7 1=1,4 
SUM=0. 

DO 8 J= 1 * 4 

8 SUM=SUM+P1CI ,JI*X1IJ) 

7 X ( I , K ) =SUM 
RETURN 
END 


SUBROUTINE EQ73(K) 

C SUBROUTINE EQ73 THIS SUBROUTINE CALCULATES THE SOLUTION VECTOR AT 
C THE STATIONIK) , GIVEN THE SOLUTION AT K*l. 

COMMON/BL 7/ PEE ( 4,4, 302 ) ,X(4,502) 

DIMENSION Z ( 4, 502 ) 

EQUIVALENCE ( X ( 1 , 1 ) ,Z ( 1 , 1 ) ) 

DO 1 1=1,4 
SUM=0. 

DO 2 J= 1 , 4 

2 SUM=SUM+PFE( I, J,K)*Z( J,K+1) 

1 Z( I ,K)=X( I ,K)-SUM 
RETURN 
END 
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SUBROUTINE F I NAL ( NMAX , I NU5 , 1N04 ) 

C CALCULATION OF SOLUTION VECTUR ASSOCIATEO WITH LAST STATION 

COMMON/BL3/NU,LAM,N,EALSIG, CH A R ,DE L / BL 5/H ( 4 ,4) ,FF(4) , J A Y ( 4 , 4 ) / BL7/ 
1PEE!4,4,502),X( 4 , 5 02 ) /BL 1 1/OMEGL ( 4 ,4 ) , CAP LL ( 4, 4) , ELU 4 ) 

DIMENSION IP I VOT ( 4) ,INUEX(4,2) ,A1(4,4),A2(4,4),A3(4,4),PSI(4,4), 
IGM(4,4) ,ETA(4) ,B(4,4) 

REAL J AY , N , NU , L AM 
21 F0RMAT(4E16.8) 

170 FORMAT! 1H0,49H THE BOUNDARY CONDITIONS AT S=SMAX ARE AS FOLLOWS//) 
169 FORMAT! 1P4E l 2 .4 ,6X 1P4E12 . 4 , 6X 1 PE 1 2. 4 ) 

167 FORMAT! 2 1 X5H0MEGA49X5HL AMDA34X 3HELL /// ) 

DO 8 1=1,4 
ELL! I ) =0 . 

DO 8 J= 1 , 4 
CAPLL ! I , J ) =0. 

8 OMEGL l I , J ) =0 . 

IF! I ND5 ) 9,9,10 

10 IF! IND5-2 ) 9,11,11 

11 CALL POLE ( N , DEL , PS 1 , GM ) 

WRITE16, 170) 

WRITE(6,40) 

40 FORMAT! 1H0, 52H THE CONDITIONS FOR A SHELL POLE HAVE BEFN GENERATE 
ID///) 

GO TO 16 

9 IF ( I ND4 ) 13,13,14 

13 READ (5,21 ) OMEGL ! 1 , 1 ) , OMEGL ( 2 , 2 ) , OMEGL ( 3 , 3 ) , OMEGL ( 4, 4 ) 

READ! 5,21) C APLLl 1 , 1 ) , CAPLL ( 2 , 2 ) , C APLL ( 3 , 3 ) ,CAPLL!4,4) 

READ (5,21) ELL 

GO TO 15 

14 READ! 5,21) OMEGL , CAPLL, ELL 

15 WRITE(6, 170) 

WR I T E (6,167) 

DO 171 1=1,4 

171 WRITE! 6, 169) OMEGL! I , 1) , UMEGL! I ,2) .OMEGL! I ,3 ), OMEGL (I ,4) , 

1 CAPLL! I, 1) .CAPLL! 1,2) .CAPLL! 1,3 ), CAPLL! 1,4), 

2 ELL ( I ) 

DO 1 1=1,4 

DO 1 J= 1 , 4 

A 1 ( I , J ) = J AY ( I » J ) /2 • +H ( I » J ) /DEL 
1 A2( I,J)=JAY( I,J)/2.-H(I, J)/OEL 
DO 2 1=1,4 
DO 2 J = 1 , 4 
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S2=0 • 

S3=0. 

DO 3 L= 1 » 4 

S2=OMEGL( I ,L)*A1CL»J) +S2 

3 S3=OMEGL 1 1 »L ) *A2 ( L » J ) +S3 
PS I ( I, J)=S3+CAPLL( I,J)/2. 

2 GM< I , J)=S2+CAPLL( 1, J)/2. 

16 DO 4 1=1,4 
DO 4 J=l, 4 
S1 = 0. 

DO 5 L= 1,4 

5 S1=S1+PSI ( I »L)*PEfc(L,J» NMAX-1 ) 

4 B( I, J)=GM( I, Jl— SI 
DO 6 1=1,4 

S1=0. 

S 2=0 . 

DO 7 J = 1 » 4 

S1=S1+PSI ( I NMAX-1) 

7 S2=S2+OMFGL( I , J)*FF( J) 

6 ETA(I)=ELL( I )-Sl-S2 

CALL MAT I NV( B, 4, ETA, 1 , DfcTEKM, I PIVOT, INDEX, 4 , 1 SCALE) 
DO 12 1=1,4 
12 X( I ,NMAX) =ETA( I ) 

RETURN 

END 



non 


SUBROUTINE STR ESS < FREQ , NMAX , N , IN05) 

SUBROUTINE STRESS — THIS SUBROUTINE CALCULATES THE SECONDARY QUANTITIES 
N XI, N XI THETA, Q XI, PHI, M THETA, M XI THETA, M XI THETA, AND 
N THETA AT EACH STATION ALONG THE SHELL. 

COMMON R( 50 2 ) /BL1/GAM<502 ) ,UMT( 502) ,OMXI (50 2> ,DEOMX( 502J/BL3/NU, 

ILAM,N, EAL SI G, CH AR , DEL/6L7/ PEE ( 4,4,502) ,X(4, 502 ) / BL5/H ( 4 , 4 ) ,FF(4), 
2JAY(4,4)/BL13/AK(3,4) , ALL (3,4) ,STHER(3) 

DIMENSION Z(4,502),Y(4) , DZ ( 4) 

EQUIVALENCE (Z(1,1),X(I,1>) 

INTEGER FREQ 
REAL N , NU , J AY , L AM 
KQUNT=0 
DO 9 1=1, NMAX 

IF ( ( I.EQ. l).OR.(I.EQ.NMAX) ) GO TO 1 
I F ( I — 2 ) 13,13,14 
14 KQUNT=KGUNT+ 1 

IF( KOUNT-FREQ) 9, 12, 12 

12 KOUNT=0 

13 DO 3 L= 1 , 4 
Y(L)=Z(L, I) 

3 DZ(L)=(Z(L, I+1)-Z(L,I-1) )/2./DEL 
GO TO 2 

1 I F ( I— 1)4, 4, 5 

4 K = 2 

GO TO 6 

5 K=NMAX 

6 DO 11 L = 1 , 4 

Y ( L ) =• 5* ( Z(L,K)+Z(L, K— 1 ) ) 

11 DZ(L)=(Z(L,K)-Z(L, K- 1 ) ) / DEL 

2 CALL HFJ( I , IND5,NMAX, 1. ) 

CALL KLT ( I , I ND5 , NMAX ) 

DO 7 L= 1 , 4 
SUM 1 = 0 • 

SUM 2 = 0 • 

D08M= 1 f 4 

SUM 1=SUM1 + H ( L , M ) *DZ ( M ) 

8 SUM2=SUM2+JAY(L ,M)*Y(M) 

7 PEE(L,2,I )=SUM1+SUM2+FF(L) 

DO 9 L= 1 , 3 

SUM 3 = 0 . 

SUM4=0. 

DO 10 M= 1 , 4 

SUM3=SUM3+AK (L, M)*OZ (M) 

10 SUM4=SUM4+ALL (L ,M)*Y(M) 

PE E(L, 1,1 )=SUM3+SUM4+STHER( L) 

9 CONTINUE 
RETURN 
END 



SUBROUTINE KLT< K, IND5.NMAX) 

C SUBROUTINE KLT THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE MATRICES 
C WHICH ALLOW THE CALCULATION OF THE QUANTITIES M-THETA , N-THETA ,AND M-XI THETA 
C AT THE STATION SPECIFIED BY THE INDEX K. 

COMMON R( 502)/BL1/GAM<502) ,UMT(502) , OMXI ( 50 2) , OEOMXI 502) /BL3/NU, 
1LAM,N,EALSIG,CHAR,DEL/BL13/AK( 3, 4), ALL (3,4) ,STHER(3) 

REAL NUtLAMtN 

CALL BDB(K,DELtNU.BtD6,D,00) 

D1=D*( l.-NU) 

D2=D*( l.-NU**2) 

0X=0MX I ( K ) 

EAL=EALSIG 
H=HHT( K.t DEL ) 

HRB = HR AI K » DEL ) 

TEMPER=TEMP(K»OEL) 

DELTAT=DELT( K,DEL) 

STHER ( 1 ) =— CHARGE AL* DEL TAT *H**3*( l./(2.«HRB )-l./HRB**2 
l+2./(3.*HRB**3) ) 

STHERI 2)=-2.*H*EAL*TEMPER /( ( l.-NU)*HRB) 

STHER(3)=0. 

I F ( IND5 ) It 1» 2 

2 I F I ( ( IND5-2) .LE.O) .AND. (K.EQ. 1) ) GO TO 8 
IF( < ( IND5-2) .GT.O) . AND. (K.EQ.NMAX) ) GO TO 8 
1 R A=R ( K ) 

GA=GAM( K ) 

0T=0MT ( K ) 

REX=3.*0T-0X 
RXE=3.*0X— OT 
RAN=N/RA 
RAN2=RAN**2 
AK ( 1 1 1 ) =0 . 

AK ( I » 2) =0 . 

AK ( 1 1 3 ) =-GA*D2 
AK ( 1 » A ) =0 . 

AK<2,1)=B*NU 
AK ( 2 » 2 ) =0 . 

AK ( 2 » 3 ) =0 . 

AK ( 2 » A ) =0 . 

AK ( 3 * 1 ) =0 . 

AK(3,2)=D1*REX/A. 

AK ( 3 t 3 ! =R AN*D 1 
AK ( 3 1 A ) =0 . 
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noon 


ALL ( 1, 1)=GA*0X*D2 

ALL ( 1, 2)=D2*RAN*OT 

ALL(1,3)=D2*RAN2 

ALL ( 1 » 4) = NU 

ALU 2, 1 ) = B*GA 

ALL ( 2 » 2 ) = B*R AN 

ALU 2,3)=B*(QT + NU*OX) 

ALL(2,4)=0. 

ALL (3, l)=-Dl*RAN*RXE/4. 
ALL(3,2)=-GA*Dl*REX/4. 

ALL ( 3 , 3 ) =— GA*RAN*D1 
ALL ( 3 » 4 ) =0. 

GO TO 6 
8 DO 3 1=1,3 
00 3 J= 1 » 4 
AK { I , J ) =0 . 

3 ALL ( I , J ) =0 . 

Cl=(N**2/2.-l.) /( l.+NU-N**2*NU/2. ) 
AK { 1, 1)=D2*0X*< (l.+NU)*Cl«-l. J 
AK( 1,2>=D2*0X*N*(NU*C1+1.) 

ALL ( 1 , 4 ) = ( NU— ( l.-NU**2)*Cl) 

AK { 2 » 1 ) =B* ( l.+NU) 

AK ( 2 * 2 ) =B*N 

ALL(2,3)=B*I l.+NU) *UX 

C2=2.+2.*NU-NU*N**2 

AK ( 3 , 1 ) =0 1*N* ( 1./C2-UX/2.) 

ALL ( 3, 4) =— N*( l.-Nu)/C2 
6 RETURN 
END 


SUBROUTINE BOB ( K , UtL , NU , B , OB , 0 ,DD > 

SUBROUTINE BDB — THIS SUBROUTINE CALCULATES THE BENDING STIFFNESS 
,D, THE MEMBRANE STIFFNESS, 6, AND THE DERIVATIVES OF D AND B , DD AND 
DB, RESPECTIVELY, FOR A SHELL COMPOSED OF A CORE HAVING NO STIFFNESS 
AND TWO SYMMETRICAL CUVtR PLATES 
REAL NU , N , L AM 
HRB=HRA(K,DEL> 

DHR8=DHRA I K, DEL > 

H=HHT ( K ,DEL ) 

DH=DHHT (K,DEL) 

D2= 1 .-NU**2 
B=2»*H/(D2*HRB) 

0=H**3* ( 3 . / I 2.*HRB)-3./HRB**2+2./HRB**3) /(D2*3. ) 

DB=2.*DH/(D2*HRB)— B*DHRB/HRB 

DD=3.*DH*D/H+H**3*UHRB/(D2*HRB**2)*(-.5+2./HRB-2./HRB**2) 

DD=DD/3. 

RETURN 

END 
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SUBROUTINE POLE < N , DEL > A 1 , A2 ) 


C SUBROUTINE POLE- THIS SUBROUTINE 
C A CLOSED SHELL 

DIMENSION A 1 ( 4 » A ) » A2 (4,4) 

REAL N 
DO 1 1=1,4 
DO 1 J= I » 4 
Aid ,JI=0. 

1 A2(I,J)=0. 

IFIN.EQ.O.) GO TO 2 

IF! (N.EQ. 1.) .OR.IN.EQ.-l.) > GO 

All 1,1) = . 5 

Al(2,2)=.5 

A2(l,l)=.5 

A2( 2 » 2) =• 5 

IFIN.NE.2.) GO TO 4 

A2 ( 3 , 3 ) = I . /DEL 

A1(3,3)=-1./DEL 

A I ( 4 , 4 ) =-l ./DEL 

A2(4,4)=1./DEL 

RETURN 

4 Al(4,4)=.5 
A2 ( 4 , 4 ) =. 5 
A1 ( 3, 3 ) =. 5 
A2 { 3 , 3 ) = . 5 
RETURN 

2 Al(l,l)=.5 
A1 ( 2 , 2 )=. 5 
A1(3,3)=-1./DEL 
Al(4,4)=-I. /DEL 
A21 1 , 1 ) =• 5 
A2(2,2)=.5 
A2(3,3)=I./DEL 
A2(4,4)=1./DEL 
RETURN 

3 A1 ( 2, 1 ) =• 5 
Al(2,2)=.5 
A1(1,I)=“1. /DEL 
Al(3,3)=.5 
Al(4,4)=.5 
A2(2,I)=.5 
A2(2,2)=»5 

A2 (.1 , 1 ) =1 ./DEL 

A2 C 3 , 3 ) = . 5 

A2(4,4)=.5 

RETURN 

END 


CALCULATES THE FINITENESS CONDITIONS FOR 


TO 3 
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SUBROUTINE OUTPUT ( FREQ , NMAX , I N01 » OEL , NO ) 

C SUBROUTINE OUTPUT THIS SUBROUTINE CONTROLS PROGRAM PRINTING ANO 

C PUNCHING. GEOMETRIC DATA IS PRINTED IF INOI IS NOT EQUAL TO ZERO. ANY 
C OR ALL OF THE ELEVEN OUTPUT QUANTITIES CAN BE PUNCHED BY SUITABLE 
C SPECIFICATION OF THE FIELDS OF THE NOL CARD 

COMMON R( 502 ) /BL1/GAM! 502 ) ,0MT!502> ,OMXI(50 2) ,DEOMX( 502 ) /BL7/PEEI 
14,4,502) ,X!4, 502) 

DIMENSION NOLI 11) , JJ( 11) , KK 1 1 1 ) , ESS I 502 ) ,Y0RD(502) 

EQUIVALENCE I X ( 1 , 1 ) , ESS I 1 ) ) , ( XI 1 , 127) .YORDIl) ) 

INTEGER FREQ 

IF ( IND1.NE.0) GO TO 51 

WRITE ( 6, 1 1 ) 

11 FORMAT I 1H1 , 16H R/RB 16H Z/RB 16H S/RB 

1 16H OMFGA THETA 16H OMEGA XI 16H DEOMEGA XI 16H 

2 GAMMA ///) 

12 FORMAT! 1P7E16. 8) 

ZED=0. 

S=0. 

DO 8 1=1, NMAX 
IF! 1-1)8, 8,9 
9 DEM=DEL 

IF < C I .EQ. 2) .OR. I I .EQ.NMAX) )DEM=.5*DEL 
S=S+OEM 

ARGU=DEM**2-!Rl I)-R( I— 1 > )**2 
IF C ARGU.LE.O. ) GO TO 8 
ZED=SQRT ( ARGU ) +ZED 

8 WRITE! 6, 12)R< I ) ,ZED,S,UMT! 1 ) ,OMXI < I ) ,DEOMX! I) ,GAM( I ) 

51 WRITE 16,2) 

2 FORMAT! 1H1, 1X5H S 1IH N XI 11H N THETA 11H N XT 
111H Q XI IIH M XI 1 1H M THETA 1 1H M XT 11H U X 
21 11H U THETA IIH W 11H PHI XI ///) 

7 FORMAT! 1XF6.3,1P11E11.4) 

K0UNT=0 

DO 33 1=1, NMAX 

IF! I I .EQ. 1) .OR. ! I .EQ.NMAX) ) GO TO 103 
I F ! I — 2 ) 14,14,15 

15 KQUNT =KOUNT + 1 

IF ( KOUNT— FREQ ) 33,16,16 

16 K0UNT=0 

14 PfcE«4,3, I ) =X ! 4 , I ) 

PEE! 1,3,1 ) = X ! 1,1 ) 

PEE 12, 3, 1 )=X!2,I) 

PEE!3,3t I ) = X I 3 , I ) 

GO TO 33 

103 IF(I-l) 104,104,105 

104 K=2 
J = 1 

GO TO 106 

105 K=NMAX 
J=NMAX 

106 PEE! 4, 3, J )=.5*< X { 4 ,K ) +X ! 4 , K-l ) ) 

PEE! 1, 3, J )=.5*( X! 1,K)+X1 1,K-1) ) 

PEE! 2, 3, J )=.5*I X!2,K)+XI2,K-1> ) 

PEE! 3,3,J )=.5*!X!3,K ) +X l 3, K— 1 ) ) 

33 CONTINUE 

NOPT S= ! NMAX -2 ) /FRfcQ+2 
SMAX=DEL*FLOAT ! NMAX— 2 ) 

WRITE!6,53)SMAX 
53 FORMAT ! 9H SMAX/R8= IPE 16. 7 ) 

ESS! 1 ) =0 . 

DO 20 1=2, NMAX 
DEM=DEL 

IF! < I. EQ. 2). OR.! I. EQ.NMAX) ) DEM=DEL/2. 



20 ESS ( I ) = ESS( I - 1 ) + D t M /SMAX 
K0UNT=0 
00 3 1 = 1, NMAX 

IF( ( I .EQ. 1) .OR. ( I .EQ.2) .OR. ( I .EQ.NMAX) ) GO TO 17 
KGUNT=K0UNT+1 
IF(KGUNT-FREQ) 3,18,18 
18 K0Ul\ir=0 

17 WRITE{6,7)ES$(I),PE£(i,2,I),PEE(2,l,I) »PEE( 2,2,1 ),PEE(3,2,I) , 

1PEE(4,3,I ),PEE< 1,1,1 ), PEE (3,1, I), PEE i 1,3,1 ) ,PE E < 2 , 3 , I > , PE E ( 3 , 3 , I ) , 
2PEE (4,2,1 ) 

3 CONTINUE 

READ (5,29) < NOL( 1 ) , 1=1 , 11 ) 

29 FORMAT (1114) 

M=0 

DO 30 1=1,11 

30 M=NOL { I ) + M 
IF(M.FQ.O) GO TO 13 
JJ ( 1) =1 

JJ( 2) =2 
J J ( 3 ) =2 
J J ( 4 ) =3 
J J ( 5 ) =4 
J J ( 6 ) =1 
JJ< 7 ) =3 
J J ( 8 ) =1 
J J ( 9 ) =2 
JJ( 10 ) = 3 
JJ< 1 1 ) =4 
KK ( 1 ) =2 
KK(2) =1 
KKI3) =2 
KK(4) =2 
KK ( 5 ) =3 
KK ( 6 ) =1 
KK ( 7 ) =1 
KK { 8 ) =3 
KK ( 9 ) =3 
KK ( 10 ) = 3 
KK ( 1 1 ) = 2 
DO 50 L = 1 ,11 

IF (NOL(L).EQ.l) GO TO 31 
GO TO 50 

31 NOL { L ) =0 
J= J J I L ) 

K=KK( L ) 

KOUNT = 0 

DO 19 1=1, NMAX 
M=I 

IF( ( I .EQ. 1) .OR. ( I .EQ.2) ) GO TO 19 
M=(NMAX-2)/FREQ+2 
I F ( I. EQ.NMAX) GO TO 19 
K0UNT=K0UNT+1 
I F ( KOUNT-FREQ ) 19, 21,21 
21 KQUN T = 0 

M=( I -2 ) /F REQ+ 2 
19 YORO(M)=PEE( J,K, I) 

Nl=( NMAX-2) /FREQ+2 

101 FORMAT (5E12.5,4X2I4) 

DO 102 1=1, Nl, 5 

102 PUNCH 101 , YORD ( I ) , YORDt I+1),Y0RD< 1 + 2) ,YQRD( 1+3) ,YORD< 1+4) ,NO,L 
50 CONTINUE 

13 RETURN 
END 
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MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS 

SUBROUTINE M AT I NV C A , N ,B , M , OE TE RM, I P I VOT , I ND EX , NMAX , I SCAL E ) 

01 MENS I UN IP I VOT (N) «A(NMAX>N)« 8<NMAX,M) , I ND EX < NM AX , 2 ) 
EQUIVALENCE (IROW,JROW), < ICUL UM , JCOLUM ) , (AMAX, T, SWAP) 

INITIALIZATION 

5 I SC AL E = 0 

6 Rl=10. 0**18 

7 R2=i.0/Rl 
10 DETERM= 1,0 
15 00 20 J=1,N 
20 I P I VOT ( J ) =0 
30 DO 550 1 = 1, N 

SEARCH FOR PIVOT ELEMENT 

40 AMAX = 0 • 0 


45 

DO 

105 J=1,N 




50 

IF 

( I P 1 VOT ( J ) - 

-1) 

60, 

105, 6 C 

60 

DO 

100 K=1,N 




70 

IF 

( I P I VOT ( K ) - 

-1) 

80 , 

100, 740 

80 

IF 

( ASS(AMAX)- 

- ABS ( A l J 

, K ) ) ) 85 , 


85 I ROW= J 
90 ICQLUM=K 
95 AM AX = A ( J , K ) 

100 CONTINUE 

105 CONTINUE 

IF (AMAX) 110,106,110 

106 DETE RM = 0 • 0 
I SC AL E = 0 
GO TO 740 

110 IP I VOT ( IC0LUM) = IPIVGT ( I COLUM ) + 1 

INTERCHANGE ROWS TU PUT PIVOT ELEMENT ON DIAGONAL 

130 IF ( IROW-ICOLUM) 140, 260, 140 

140 DETERM=-QETERM 
150 DO 200 L= 1 , N 
160 SW AP = A ( I ROW , L ) 

170 A( IRQW,L)=A (I COLUM , L ) 

200 A( ICQLUM,L)=SWAP 
205 I F ( M ) 260, 260, 210 
210 DO 250 L= 1 , M 
220 SW AP = B ( I ROW , L ) 

230 B( IROW,L)=B( I COLUM , L ) 

250 B ( I COLUM ,L)=SWAP 
260 I NOEX ( l , 1 ) = I ROW 
270 INDEX(I,2)= I COLUM 
310 PIVOT=A( ICQLUM, ICULUM) 

SCALE THE DETERMINANT 

1000 PIVOT I=PI VOT 

1005 IFC ABS(DE TERM) -Rl) 1030,1010, 1010 
1010 DETERM=UETERM/R1 
ISCALE=ISCALE+1 

IF (ABS(DETFRM)-Rl ) 1060,1020,1020 
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1020 D£TERM=0ETERM/R1 
ISCALE=ISCALE+1 
GO TO 1060 

1030 IF tABStDE TERM )-R2> 1040,1040,1060 
1040 0ETERM=DETERM*R1 
ISCALE=ISCALE— 1 

IF (ABSl 0ETERM1-R2) 1050, 1050,1060 
1050 DETERM=DETERM*R1 
ISCALE=ISCALE— 1 

1060 IF (A8 St PIVOT I )-Rl) 10 90,1070,1070 
1070 PIVOTI=PI V0TI/R1 
ISCALE=ISCALE+1 

IFtABStPI VOTIl-Rl 1320,1080, 1080 
1080 PIVOTI=PI V0TI/R1 
ISCALE=ISCALE+1 
GO TO 320 

1090 IFtABStPI VOTI )-R2 >2000,2000, 320 
2000 PIVUTI=PIV0TI*R1 
ISCALE=ISCALE— 1 

IFtABStPI V0TIJ-R2) 2010, 2010, 320 
2010 PIV0TI=PIV0TI*R1 
ISCALE=ISCALE— 1 
320 DETERM=DE TERM*PIVOT I 
C 

C DIVIDE PIVOT ROW BY PIVOT ELEMENT 

C 

330 At ICOLUM, ICOLUM)=1.0 
340 DO 350 1=1, N 

350 At ICOLUM, L)=A( ICOLUM.L) /PIVOT 
355 IF(M) 380, 380, 360 
360 DO 370 L=1,M 

370 B( ICOLUM, L ) = B t ICOLUM, L> /PI VUT 
C 

C REDUCE NON-PIVOT ROWS 

C 

380 DO 550 L 1 = 1 , N 
. 390 IF l L 1- 1 COLUM > 400, 550, 400 
400 T=A(L1, ICOLUM) 

420 AtLl, IC0LUM)=0.0 
430 DO 450 L= 1 , N 

450 A(Ll,L)=AtLl,L>— A(I CULUM , L ) *T 
455 IF(M) 550, 550, 460 
460 DO 500 L= 1 , M 

500 Bt LI ,L)=B(L1,L)-Bt ICOLUM, L)*T 
550 CONTINUE 
C 

C INTERCHANGE COLUMNS 

C 

600 DO 710 1=1, N 
610 L=N+ 1— I 

620 IF t INDEXtL,l)-INDEXlL,2>) 630, 710, 630 
630 JROW=INDEXtL,l) 

640 JC0LUM=INDEXIL,2> 

650 DO 705 K= 1 , N 
660 SWAP=AtK, JROW) 

670 At K, JRGW) = A(K , J COLUM ) 

700 At K, JCOLUM )=SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 
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APPENDIX C 


MEMBRANE AND BENDING STIFFENESSES 

The membrane and bending stiffnesses of the shell are defined by equations (A24) 
and (A25), respectively. The membrane and bending stiffnesses of a sandwich shell such 
as that shown in sectional view by figure 5 can be used in the program. For such a 
sandwich, shear deformations in the core are neglected as well as core compressibility 
through the thickness, yet the core is assumed to have no bending or membrane stiff- 
ness. The face sheets are equal and assumed to carry only membrane forces. These 
assumptions lead to the expression: 

2Et 


b = 


Eohol 1 - v 2 ) 


(Cl) 


Let 


then equation (Cl) becomes 


similarly, 


h 

V = r- 


6 -r 


or 


b = ?e a 

e 0 (i - t ) 5 


d " h2t ■ 3 “ 2 + 2 t 1 ft 


(C2) 


(C3) 


„ 1 „3 /3 3 . 2 \ E 

d = — V o* " -O + “9 ■ 


3(l - t) \ 25 5 2 63^0 


(C4) 


The derivatives of b and d are also required and can be obtained from equations (C3) 
and (C4) 


b’ = , 2E 0 , (v' - 5 5') 

Eo(i - t * 2 )6 ' 5 > 


(C5) 


d' = ^ 77 'd + 
V 


77 3 5’ 


1 2 2 \ E 

(1 - v 2 )& 2 \ 2 + 6 6 2 / e o 


(C6) 
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APPENDIX C 


Note that for the isotropic case b and d reduce to the well-known relations 

b 1 h E 

(l - „2) h Q E 0 


1 E fh 


12(l - i/2) E o \h Q/ 
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APPENDIX D 


THERMAL FORCES AND MOMENTS 


Using the definition of thermal force and moment given by equations (A26) and 
(A27) with the assumed form of the temperature distribution equation (31) gives the fol- 
lowing equations: 


t 


(n) _ 


2EoT^ 




1 

u) 6 


(Dl) 



aEaAT 

3"„U ' 


Lzzi 

v) 8 


I- 



Differentiating equations (Dl) and (D2) gives 


(D2) 


t (n)’ _ 

l T “ 


2Ea 

CT 0 (i - u) 6 




VAn) 
5 T 


(D3) 


m 


(n)’_ 
T ~ 


aEa 

3a 0 (l - v ) 




ATjV 


+ 3?7^r]'ATi] + ATj 



— + — 
2 6 



(D4) 
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APPENDIX E 


OUTPUT LISTING FOR CLOSED SPHERICAL SHELL SEGMENT 

The following is the output list for the problem of a spherical segment subjected to 
a hydrostatic pressure and having a horizontal elastic boundary restraint. The problem 
is described in the text in the section entitled "Example Problems." 
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PROBLEM NUMBER 


1 


SPHERICAL CAP WITH HORIZONTAL EDGE RESTRAINT SUBJECT TO HYDRO PRESSURE 

INPUT DATA 


1 •OOOOOOOE 03 
3.0000000E-01 
0.0000000E-39 
I. OOOOOOOE 00 
0.0000000E-39 
52 


REFERENCE LENGTH = 

POISSONS RATIO 
TEMP COEFFICIENT = 

REF. THICKNESS 
FOURIER INDEX 
NUMBER OF STATIONS^ 

THE INDICATORS ARE SET AS FOLLOWS 
I ND1= 0 IND2= 0 I ND3= 


1 I ND4- 1 IND5= I 


THE BOUNDARY CONDITIONS AT S=0 ARE AS FOLLOWS 

THE CONDITIONS FOR A SHELL POLE HAVE BEEN GENERATED 


THE BOUNDARY CONDITIONS AT S=SMAX ARE AS FOLLOWS 


OMEGA 


LAMDA 


ELL 


-O.OOOOE-39 
-0.0000E-39 
8 « 6603E-01 
-0.0000E-39 


-0# 0000E-39 
-O.OOOOE-39 
-O.OOOOE-39 
-O.OOOOE-39 


-0. 0000E-3 9 
— 0. OOOOE-39 
5.0000E-01 
-0.0000E-39 


-0. OOOOE— 39 
-O.OOOOE-39 
— 0.0000E-39 
-O.OOOOE-39 


-5.0000E-0I 
0.0000E-39 
8 • 6603E-02 
-O.OOOOE-39 


-O.OOOOE-39 
l.OOOOE 00 
-0.0000E-39 
-O.OOOOE-39 


8.6603E-01 

O.OOOOE-39 

5.0000E-02 

-O.OOOOE-39 


-0.00006-39 
-0.0000E-39 
-O.OOOOE-39 
l.OOOOE 00 


-0.0000E-39 

-0.0000E-39 

-O.OOOOE-39 

-O.OOOOE-39 


05 

CO 


APPENDIX E 



05 


R/RB 


Z/RB 


S/RB 


OMEGA THETA OMEGA XI 


DEOMEGA XI 


GAMMA 


0* 000000 OOE— 39 
5 • 23596370E-03 
1.57073170E-02 
2 • 61 769480 E- 02 
3 • 6643 70 8 OE— 02 
4.71064500E-02 
5.75640270E-02 
6.80152890E-02 
7 • 845909 50E-0 2 
8.88942960E-02 
9.93197490E-02 
1 • 097343 1 0E-01 
1.20136 840E-01 
!• 305261 90E-01 
1 .40 90 1230 E-0 1 
1.51260820E-01 
1 . 6160 3 8 20 E-0 1 
1 • 71929 100E-0 1 
1.82235520E-01 
1.92521970E-01 
2. 02 787 2 90 E-0 1 
2.1 3030380E-01 
2.2325011 0E-01 
2 . 33445 360E-0 1 
2.4361501 OE— 0 1 
2 • 53 757940E-0 1 
2.63873050E-01 
2. 73959220 E- 01 
2. 84015340E-01 
2 .94040 3 20 E-0 1 
3 • 04033 060 E-0 1 

3. 13992450E-01 
3 • 2391 7410 E-0 1 
3.33 806 85 OE- 01 
3.43659690E-01 
3 . 53474840E-01 
3. 6325 1 230E-0 1 
3. 729 87 780 E- 01 
3. 82 683430 E-0 1 
3.92337110E-01 
4. 01947770E-01 
4.11 514350E-01 
4. 21 0358 10 E-0 1 
4.3051 1090E— 0 1 
4 .39939 160 E-0 1 
4.49318990E-01 
4. 586495 5 0E-01 
4* 679298 10E -01 
4 • 77158760E-01 
4.86335370E-01 
4. 95458660E— 0 1 
4. 99 999990 E-0 1 


0.000000 00 £—39 
1 • 5 85 080 20 E-0 5 
1 . 29997290 E— 04 
3.5 1 582440 E— 04 
6.82035750E-04 
1.12 168680E— 03 
1 . 6 70657 10 E— 03 
2.32 896940E-03 
3.09655760E-03 
3.97340650E-03 
4.9 5 940 750 E— 03 
6 . 05447 100 E -03 
7. 258508 90 E-0 3 
8.57137240E-03 
9. 99 2 92 8 40 E- 03 
1. 15230340E— 02 
1.31 6 15020E-02 
1.49081790E-02 
1. 676286 60E-02 
1 • 87253640E— 02 
2.07954630E-02 
2 . 2 9729340 E— 02 
2. 525753 00 E -02 
2 . 76490 1 30 E -02 
3.01 47 12 00 E— 02 
3.275157 10E— 02 
3.54620880E-02 
3. 82783790 E-02 

4. 12001 2 10E— 02 
4 .42 270220 E-02 
4. 73 5 869 80 E— 02 
5.05948640E— 02 
5 * 3935 1300E-02 
5 .737914 90 E-02 
6. 092652 30 E— 02 
6.457689 10E— 02 
6. 83 29 8 1 60 E— 02 
7.21 849260E— 02 
7. 614176 50 E— 02 
8.01999340E— 02 
8.43589570E-02 
8.861 83 8 80 E— 02 
9. 297776 50 E— 02 
9. 74366070E— 02 
1 • 0 1994430 E-0 1 
1 .06 6 5 07 20 E-0 1 
1. 1 1404980E-0 1 
1.162566 90 E-01 
1.2120 53 00 E-01 
1 • 26 2 503 00 E— 01 
1.31391110E-01 
1.33997220E-01 


O.OOOOOOOOE-39 
5 • 23598770E-03 
1 • 57079 63 OE— 02 
2 * 61 799390E— 02 
3. 66 5191 40 E— 02 
4 .71238890 E-02 
5 • 7595865 OE -02 
6. 80678400E— 02 
7 • 85398160E-02 
8*901 17910E— 02 
9 . 94837670E-02 
1 « 09955740E-01 
1.20427720E-01 
1 • 30899690E— 01 
1.41371670E-01 
1 • 51 843640E— 01 
1 . 62315620E-01 
1 . 72787590E-01 
1. 83259570E-01 
1 . 93731 540E -01 
2 • 042 03 5 2 OE -01 
2 • 14675500E-01 
2 • 25 147470E-0 1 
2.35619450E-01 
2 .46091420E— 01 
2.56563400E-01 
2 • 67035370E— 01 
2 • 77507340E— 01 
2 . 87979320E-01 
2 • 9845 1290E-01 
3 • 08923260E-01 
3. 19395240E— 01 
3 • 298672 10E-01 
3 . 40339 1 80E-01 
3.50811160E-01 
3.61283130E— 01 
3.717551 OOE-Ol 
3 • 82227080E-01 
3. 92699050E-01 
4 • 031 71030E-01 
4. 13643000E— 01 
4. 24114970E— 01 
4. 34586950E— 01 
4.45058920E-01 
4.55530890E-01 
4.66002870E-01 
4. 76474840E— 01 
4 • 86946810E— 01 
4 • 9741 8790E— 0 1 
5.07890760E-01 

5. 18 362730E-01 
5 • 23598720E-01 


1 * OOOOOOOOE 00 
1.0 OOOOOOOE 00 
1. OOOOOOOOE 00 
l. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
l. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
l. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 


1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
i. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1, OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 
1. OOOOOOOOE 00 


0. OOOOOOOOE -39 

O.OOOOOOOOE-39 

0. OOOOOOOOE- 39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

0. OOOOOOOOE -39 

Q.OOOOOOOOE-35 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 

O.OOOOOOOOE-39 


O.OOOOOOOOE-39 
1.90984190E 02 
6.36567420E 01 
3 . 81 884600E 01 
2 • 72714870E 01 
2.12049490E 01 
1 • 7343 1540E 01 
1.46685290E 01 
1.27062050E 01 
1. 12047800E 01 
1 • 00187080E 01 
9.05788660E 00 
8 « 26355480E 00 
7 • 59575420E 00 
7.026366306 00 
6.53502940E 00 
6. 10663600E 00 
5.72974170E 00 
5. 39551720E 00 
5.09704260E 00 
4. 828817406 00 
4. 58641420E 00 
4. 36622930E 00 
4.16529980E 00 
3.981 16690E 00 
3.81177330E 00 
3.65538440E 00 
3. 51052730E 00 
3.37594350E 00 
3. 25055090E 00 
3.13341410E 00 
3.02372070E 00 
2.92076 1 OOE 00 
2.82391290E 00 
2.73262840E 00 
2.64642330E 00 
2.56486740E 00 
2.48757820E 00 
2.41421360E 00 
2 • 34446720E 00 
2.27806360E 00 
2 • 21475450E 00 
2. 15431560E 00 
2. 09654360E 00 
2. 04125400E 00 
1*98827870E 00 
1.93 746450 E 00 
1 • 88867140E 00 
1.84177090E 00 
1.79664540E 00 
1.75318670E 00 
1.73205080E 00 



S N XI N THETA N XT 0 XI H XI H THETA M XT U XI U THETA W PHI XI 


SMAX/RB= 5.235987TE-01 

0. 000-5. 0000E-01-5.0000E-0I O.OOOOE-39-O.OOOOE-39 1 .9635E-02 1.9635E-02 O.OOOOE-39 O.OOOOE-39-O.OOOOE-39-5.2879E 00 O.OOOOE-39 
0. 010-5. 0011E-01-5.0008E-01-O.OOOOE-39-3.6883E-07 1.9635E-02 1 .9696E-02-0 .0000E-39 2 . 5855E- 02-0. OOOOE-39-5 . 2879E 00 8.6744E-04 
0. 030-5. 0003E-01-5.000IE-01-0.0000E-39-8.5912E-07 1.2155E-02 1 . 5403E-02-0. OOOOE-39 7.7562E-02 0. 0000E-39-5. 2873E 00 2.2163E-03 

0. 050-5. O0O1E-01-5.OOOOE-01-O.OO0OE-39-6.7961E-O7 5.9719E-03 1 .O938E-O2-0.O000E-39 1.2926E-0I 0. OOOOE-39-5. 2863E 00 2.8742E-03 

0. 070-5. OOOOE-OI- 5 .OOOOE-OI-O. 0000E-39-4. 5585E-07 1.8935E-03 7.4838E-03-O.OOOOE-39 1 .8095E-01-0. OOOOE-39-5 .2847E 00 3.0431E-03 
0. 090-4. 9999E-01-5.0000E-0 1-0. 000 0E- 39-2. 7285E-07-3 .8242E-04 5 . 059 IE-03-0 .OOOOE-39 2. 326 IE-01-0. OOOOE-39-5. 2825E 00 2.9279E-03 
0. 1 10-4. 9999E-01-5.0001E-01-0. 000 0£- 39-1. 4675E-07-1 .4043E-03 3. 45 13E-03-0. OOOOE-39 2. 8425E-01-0. OOOOE-39-5. 2799E 00 2.6795E-03 
0. 130-4. 9998E-01-5.0001E-01-O.OOOOE-39-6.8667E-08-1 .6923E-03 2. 4263E-03-0. OOOOE-39 3. 3586E-01-0. OOOOE-39-5. 2767E 00 2.4002E-03 
0. 150-4. 9998E-01-5.0002E-01 O.OOOOE-39-2.2707E-08-1 .5771E-03 1.7912E-03 O.OOOOE-39 3. 8744E-01-0. 00006-39-5 .2729E 00 2.1385E-03 

0. 170-4. 9998E-01-5.0002E-01 O.OOOOE-39 1 .2151E-09- 1 .271 5E-03 1.4140E-03 0.0000E-39 4. 3897E-01-0. OOOOE-39-5. 2686E 00 1.9228E-03 

0. 190-4. 9998E-01-5.0002E-01 O.OOOOE-39 7.7308E-09-9.2149E-04 1.1872E-03 O.OOOOE-39 4. 9045E-01-0. OOOOE-39-5. 2637E 00 1.7531E-03 

0. 210-4. 9998E-01-5.0002E-01 0.0000E-39 4 .9003E-09-6 .671 0E-04 1.0277E-03 0.0000E-39 5. 4188E-01-0. OOOOE-39-5. 2583E 00 1.6266E-03 

.0. 230-4. 9997E-01-5.0003E-01 0.0000 E-39-5.3300E-09-4.9734E-04 9.1359E-04 0.0000E-39 5. 9325E-01-0. OOOOE-39-5. 2524E 00 1.5434E-03 

0. 250-4. 9997E-01-5.0003E-0 1-0. OOOOE-39- 1. 9341 E-08-5 .3453E-04 7.678 IE-04-0 .0000E-39 6. 4456E-01-0. OOOOE-39-5. 2460E 00 1.4664E-03 
0. 270-4. 9997E-01-5.0003E-0 1-0. 0000E- 39-2. 572 2E-03-6 .9523E-04 5. 9698E-04-0. OOOOE-39 6. 9580E-01-0. OOOOE-39-5. 2390E 00 1.3758E-03 
0. 290-4. 9997E-01-5.OOO3E-O1-0.OO0OE- 39-2. 3 194E-08-S.8310E-04 4. 276 IE-04-0. OOOOE-39 7. 4695E-01-0. OOOOE-39-5. 2314E 00 1.2717E-03 
0. 310-4. 9997E-0 1-5. 0003E-0 1-0. 0000 E- 39-1. 08 25E-08-1. 001 6E-03 2. B829E-04-0. OOOOE-39 7. 9803E-01-0. OOOOE-39-5. 2233E 00 1.1570E-03 

0. 330-4. 9997E-01-5.0003E-01 O.OOOOE-39 1 .0299E-08-9 .4485E-04 2.1598E-04 O.OOOOE-39 8. 4902E-01-0. OOOOE-39-5. 2147E 00 1.0460E-03 

0. 350-4. 9996E-01-5.0003E-01 O.OOOOE-39 3 .6780E-08-6 .4659E-04 2.3962E-04 O.OOOOE-39 8. 9992E-01-0. OOOOE-39-5. 2056E 00 9.6437E-04 

0. 370-4. 9996E-01-5.0003E-01 O.OOOOE-39 8.6276E-08-7.4387E-05 3.7242E-04 0.0000E-39 9. 5072E-01-0. OOOOE-39-5. 1959E 00 9.2933E-04 

0. 390-4. 9996E-01-5.0003E-01 O.OOOOE-39 1.6799E-07 1.2081E-03 7.6398E-04 O.OOOOE-39 1.0014E 00-0. OOOOE-39-5 . 1857E 00 9.9790E-04 

0.410-4. 9995E-01-5.0003E-01 O.OOOOE-39 2.7539E-07 3.3991E-03 1.5018E-03 0.0000E-39 1.0520E 00-0. OOOOE-39-5 . 1749E 00 1.2613E-03 

0. 430-4. 9995E-01-5.0003E-01 0.0000E-39 4.0124E-07 6.7935E-03 2.7062E-03 O.OOOOE-39 1.1025E 00-0. OOOOE-39-5. 1637E 00 1.8364E-03 

0.450-4. 9994E-01-5.0005E-01 0.0000E-39 5.0389E-07 1 .1429E-02 4.4244E-03 0.0000E-39 1.1528E 00-0. OOOOE-39-5. 1519E 00 2.8686E-03 

0. 470-4. 9994E-01-5.0007E-01 0.0000E-39 4.8082E-07 1 .6736E-02 6.4969E-03 O.OOOOE-39 1.2030E 00-0. 0000E-39-5. 1396E 00 4.4493E-03 

0. 490-4. 9993E-01-5.0011E-01 O.OOOOE-39 1.4090E-07 2 .0645E-02 8.2689E-03 0.0000E-39 1.2531E 00-0. OOOOE-39-5. 1268E 00 6.5333E-03 

0. 510-4. 9994E-O1-5.O018E-O1-O.OO00E-39-7.7C39E-O7 1 .8699E-02 8. 2632E-03-0. OOOOE-39 1.3031E 00-0. OOOOE-39-5 . 1 135E 00 8.7113E-03 

0. 530-4. 9994E-01-5.O027E-O1-O.OO00E-39-2.4711E-O6 3.7115E-03 4. 0103E-03-0. OOOOE-39 1.3529E 00-0. OOOOE-39-5. 0997E 00 9.9024E-03 

0. 550-4. 9996E-01-5.0036E-01-O.COODE-39-5.1343E-06-3.3035E-02-7.6211E-03-0. 0000E-39 1.4026E 00-0. OOOOE-39-5. 0854E 00 8.1375E-03 
0. 570-4. 9999E-01-5. 004 IE-01-0. 0000 E-39- 8. 7291E-06-1 .0202E-0 1-3. 0505E-02-0. 0000E-39 1.4521E 00-0. 0000E-39-5. 0705E 00 3.7716E-04 
0. 590-5. 0004E-01-5.0035E-01-O.OOOOE- 39-1. 2598E-05-2 .1099E-0 1-6. 7849E-02-0. OOOOE-39 1.5014E 00-0. OOOOE-39-5. 0549E 00-1 .7435E-02 
0. 610-5. 0010E-01-5. 0004E-0 1-0. OOOOE-39- 1 .5010E-05-3 .5648E-01- 1. 1943E-01-0. OOOOE-39 1.5506E 00-0. OOOOE-39-5. 0387E 00-4.9564E-02 
0. 630-5. 0014E-0 1-4. 9934E-0 1-0. OOOOE-39- 1. 265 1 E-05-5 .1033E-01-1 .7706E-0 1-0. OOOOE-39 1.5996E 00-0. OOOOE-39-5. 021 4E 00-9.8425E-0.2 
0. 650-5. 0012E-01-4.9805E-01 0. 0000E-39-2 . 1 785E-07-6.0105E-01-2 . 1815E-01 O.OOOOE-39 1.6485E 00-0. OOOOE-39-5. 0031E 00-1 .6078E-01 
0. 670-4. 9998E-01-4.9607E-01 O.OOOOE-39 2.9503E-05-4.9225E-01-1.9814E-01 O.OOOOE-39 1.6971E 00-0. 0000E-39-4.9835E 00-2.2159E-01 
0. 690-4. 9960E-01-4.9353E-01 0.0000E-39 8.4371E-05 3 .3679E-02-4.4108E-02 O.OOOOE-39 1.7455E 00-0.C000E-39-4.9629E 00-2.4582E-01 

0. 710-4. 9888E-01-4.9102E-01 O.OOOOE-39 1.6956E-04 1 .2705E 00 3.4491E-01 0.0000E-39 1.7937E 00-0. OOOOE-39-4. 941 8E 00-1.6955E-01 

0. 730-4. 9771E-01-4.8995E-01 0.0000E-39 2.8097E-04 3.5352E 00 1.0824E 00 O.OOOOE-39 1.8416E 00-0. OOOOE-39-4. 9218E 00 1.0550E-01 

0. 750-4. 9609E-01-4.9292E-01 O.OOOOE-39 3.9447E-04 7.0273E 00 2.2503E 00 0.0000E-39 1.8894E 00-0. OOOOE-39-4. 9055E 00 7.0634E-01 

0. 770-4. 9428E-01-5.0393E-01 O.OOOOE-39 4.5118E-04 1.1555E 01 3.8104E 00 O.OOOOE-39 1.9371E 00-0. OOOOE-39-4. 8970E 00 1.7597E 00 

0. 790-4. 9295E-01-5.2823E-01 O.OOOOE-39 3.4093E-04 1.6096E 01 5.4596E 00 0.0000E-39 1.9848E 00-0. OOOOE-39-4. 9018E 00 3.3223E 00 

0. 810-4. 9346E-01-5.7110E-01-O.OOOOE-39-1.0912E-04 1.8188E 01 6.4259E 00-0. OOOOE-39 2.0327E 00-0. OOOOE-39-4. 9250E 00 5.2524E 00 

0. 830-4. 9807E-01-6. 3 512E-0 1-0. 0000 E-39- 1 .1280E-03 1.3265E 01 5.2378E 00-0. OOOOE-39 2.0812E 00-0. OOOOE-39-4. 9692E 00 7.0088E 00 

0. 850-5. 0999E-01-7. 1486E-01-0. OOOOE-39-2. 9543E-03-5 .7989E 00-4. 4895E-01 -0. 0000E-39 2.1303E 00-0. OOOOE-39-5. 0282E 00 7.3877E 00 
0. 870-5. 3296E-01-7.8867E-01-0. OOOOE-39-5. 7190E-03-4.8376E 01-1.3789E 01-0. OOOOE-39 2.1803E OO-O.0OOOE-39-5.0793E 00 4.2523E 00 

0. 890-5. 7000E-01-8.0773E-01-O.OOOOE-39-9.2224E-03-1 .2410E 02-3.8159E 01-0. OOOOE-39 2.2305E 00-0. OOOOE-39-5 . 0716E 00-5.6051E 00 

0. 910-6. 2075E-01-6.8469E-01-O.0O00E-39-1.2579E-O2-2.3795E 02-7.5616E 01-0. 0000E-39 2.2796E 00-0. OOOOE-39-4. 9152E 00-2.6208E 01 

0. 930-6. 7703E-01-2.8737E-01-O.OOOOE-39-1.3741E-02-3.8097E 02-1.2394E 02-0. 0000E-39 2.3248E 00-0. 0000E-39-4. 4750E 00-6.1333E 01 

0.950-7. 1624E-01 5.5 159E-01-O.OOOOE-39-8 .9860E-03-5 .1 558E 02-1.7188E 02-0. OOOOE-39 2.3610E 00-0.0000E-39-3.5820E 00-1.1209E 02 

0. 970-6. 9351E-01 1.9923E 00 0.0000E-39 7.3936E-03-5 .5591E 02-1.9261E 02 O.OOOOE-39 2.3813E 00-0. OOOOE-39-2. 0780E 00-1.7254E 02 

f 5 0.990-5.3477E-01 4.0S93E 00 O.OOOOE-39 4. 2745E-02-3 .4706E 02-1.3671E 02 O.ODOOE-39 2.3775E 00-0. OOOOE-39 8. 1551E-02-2.2306E 02 

W 1.000-3. 5297E-01 5.3445E 00 O.OOOOE-39 6.6343E-02 7 .6294E-06-3.4973E 01 0.0000E-39 2.3601E 00 O.OOOOE-39 1.3626E 00-2.4230E 02 
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TABLE 1.- DESCRIPTION OF SUBROUTINES 


FORTRAN name 

Description 

HFJ 

Calculates the elements of the H, J, and f matrices 

FORCE 

Calculates the forcing function due to surface loads and/ or temper- 
ature distribution 

ABCG 

Calculates the elements of the A, B, C, and g matrices 

INIT 

Calculates the initial P and x matrices 

EFG 

Calculates the elements of the E, F, and G matrices 

PANDX 

Calculates the P and x matrices at each station 

FINAL 

Calculates the vector z associated with station NMAX 

EQ73 

Calculates the solution vector z at each shell station 

STRESS 

Calculates explicitly the moments, forces, and displacements at 
each shell station 

OUTPUT 

Controls listing of calculated data 

MATINV 

Performs matrix inversion and matrix multiplication 

KLT 

Calculates the elements of the matrices necessary for calculation 
of the output quantities m^, t^, and m^ 

BDB 

Calculates the membrane and bending stiffnesses at each shell 
station 

POLE 

Sets up the finiteness conditions at a pole in the shell geometry 
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TABLE 2.- GLOSSARY OF FORTRAN VARIABLES 


Variable 

Subroutine where 
variable is calculated 
or defined 

Description 

A(4,4) 

ABCG 

Matrix of coefficients Aj as defined 
by equations (20) 

B 

BDB 

Membrane stiffness 

B(4,4) 

ABCG 

Matrix of coefficients Bj as defined 
by equations (20) 

C(4,4) 

ABCG 

Matrix of coefficients Ci as defined 
by equations (20) 

CAPL1(4,4) 

SHELLS 

Matrix [a]j, see equation (15) 

CAPLL(4,4) 

SHELLS 

Matrix [a]^, see equation (15) 

CEE(4) 

FORCE 

Load vector e^ as defined by equation (13) 

CHAR 

SHELLS 

a, characteristic shell dimension 

CONST (100) 

SHELLS 

Unassigned storage 

D 

BDB 

Bending stiffness 

DB 

BDB 

Derivative of membrane stiffness 

DD 

BDB 

Derivative of bending stiffness 

DEL 

SHELLS 

A£, arc length increment 

DEOMX(502) 

INPUT 

dco^/dl 

DH 

BDB 

n' 

DHRB 

BDB 

6' 

E(4,4) 

EFG 

Matrix QeJ^ see equation (13) 

EALSIG 

SHELLS 

E aj a 0 , see equation (Dl) 

ELI (4) 

SHELLS 

Vector { l see equation (15) 

ELL(4) 

SHELLS 

Vector {Y}^, see equation (15) 

F(4,4) 

EFG 

Matrix [f]j, see equation (13) 

FF(4) 

HFJ 

Vector (f), see equation (21) 

G(4,4) 

EFG 

Matrix [gJ^, see equation (13) 


68 



TABLE 2.- GLOSSARY OF FORTRAN VARIABLES - Concluded 


Variable 

Subroutine where 
variable is calculated 
or defined 

Description 

GAM (502) 

INPUT 

r 

H 

BDB 

V 

H(4,4) 

HFJ 

Matrix [h)j, see equation (21) 

HRB 

BDB 

5 

IND1 

SHELLS 

Logical control for output of geometrical 



data 

IND2 

SHELLS 

Logical control of data into CONST 

IND3 

SHELLS 

Logical control for calculation of geometri- 



cal data 

IND4 

SHELLS 

Logical control for boundary condition 



matrices 

IND5 

SHELLS 

Logical control for closed shell 

JAY (4, 4) 

HFJ 

Matrix Qf]p see equation (21) 

LAM 

SHELLS 

X 

N 

SHELLS 

n (Fourier coefficient number) 

NMAX 

SHELLS 

N, total number of shell stations 

NU 

M 

SHELLS 

V 

OMEGl(4,4) 

SHELLS 

Matrix (fClp see equation (15) 

OMEGL(4,4) 

SHELLS 

Matrix jjf] see equation (15) 

OMT(502) 

INPUT 


OMXI(502) 

INPUT 


PEE(4,4, 502) 

INIT, PANDX 

Matrix jjpj p see equation (25) 

R(502) 

INPUT 

P 

SMAG(4) 

ABCG 

Vector defined by equation (20) 

SMAX 

OUTPUT 

Total nondimensional arc length 

X(4,502) 

INIT, PANDX 

Vector {x^p see equation (25) 
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TABLE 3.- COMMON LOCATION OF VARIABLES REQUIRED 
FOR USER- PREPARED SUBPROGRAMS 


1 

Variable name 

COMMON location 

CONST (100) 

BL14 

R(502) 

Blank COMMON 

GAM(502), OMT(502), OMXI(502), 
DEOMX(502) 

BLl a 

NU, LAM, N, EALSIG, CHAR, 

BL3 a 

DEL 



a Variables in this COMMON block must be in 


order specified. 



TABLE 4.- RELATIONSHIP BETWEEN NONDIMENSIONAL PROGRAM 
OUTPUT AND PHYSICAL QUANTITIES 


Relationship 

Storage location 
at output time 

N ondim ensional 
output variable name 


(a) 


£ a o h o 

PEE(1,2,K) 

N XI 

t (n) . N i n> 
6 CT o h o 

PEE(2,1,K) 

N THETA 

c-t- 

imo 

ii 

Q Zl 
O .m-g- 

PEE(2,2,K) 

N XT 

,(n) _ 

£ ^O h O 

PEE(3,2,K) 

Q XI 

m< n > = M<"> a , 
5 4 "oXo 

PEE(4,3,K) 

M XI 

o u o 

PEE(1,1,K) 

M THETA 

m (n) - M^ n) a 
£0 M ^e h 3 
u o n o 

PEE(3, 1,K) 

M XT 

w°ls° 

s> 

ii 

JUJ' 

3 

PEE(1,3,K) 

U XI 

(n) _ yOi) ^o 
u e " u 0 acr Q 

PEE(2,3,K) 

U THETA 

w (n) = wW ^ 2 . 

ao o 

PEE(3,3,K) 

W 

(n) ,(n) E o 

n ■ ®o 

PEE(4,2,K) 

PHI XI 


a Where K is an index specifying the shell station; 
K = 1, 2, . . NMAX. 
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TABLE 5.- OUTPUT FORMAT 


PROBLEM NUMBER (14) 

(PROBLEM DESCRIPTION) 

INPUT DATA 


REFERENCE LENGTH 

= (1PE16.7) 

POISSON' S RATIO 

= (1PE16.7) 

TEMP PARAMETER 

= (1PE16.7) 

REFERENCE THICK 

= (1PE16.7) 

FOURIER INDEX 

= (1PE16.7) 

NUMBER OF STATIONS 

= (14) 

THE INDICATORS ARE SET AS FOLLOWS 

INDl = (14) 

IND2 = (14) 


IND3 = (14) 

THE FOLLOWING CONSTANTS APPEAR IN COMMON BLOCK BL14 


(Data are printed columnwise according to FORMAT 1X1P4E16.8) 
THE BOUNDARY CONDITIONS AT S = 0 ARE AS FOLLOWS 
OMEGA LAMDA 


IND4 = (14) IND5 = (14) 


ELL 


(Matrix elements are printed by rows using FORMAT 1P12.4) 

THE BOUNDARY CONDITIONS AT S = SMAX ARE AS FOLLOWS 

OMEGA LAMDA ELL 

(Matrix elements are printed by rows using FORMAT 1P12.4) 

R/RB Z/RB S/RB OMEGA THETA OMEGA XI DEOMEGA XI GAMMA 

(1PE16.8) (1PE16.8) (1PE16.8) (1PE16.8) (1PE16.8) (1PE16.8) (1PE16.8) 

S N XI N THETA N'XT Q XI M XI M THETA M XT U XT U THETA W PHX XI 

(1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) (1PE11.4) 





f c ) Moment resultants. (d) Loads per unit area. 

Figure 2.- Positive sense of forces, moments, and loads on shell segment. 
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Figure 8.- Variation of axial force Fv with increment size for short cone (S/a = 0.0628) having nose half-angle of 7r/3 radians. Symbols 
represent results of numerical calculation; N r mav is the value of meridional stress calculated for A/(S/a) = 0.1. 






Figure 11.- Variation of with arc length for hydrostatically loaded spherical segment having horizontal elastic restraint at shell edge. 

(t> n = — — — = 0.1 where k is a proportionality constant relating horizontal edge force and edge displacement. 

0 6 h 0 E 0 
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